methacetin-fitting
Daniel Lill
23 November 2017
Load all important libraries
library(dMod)
library(stringr) # Um bequem mit strings zu arbeiten
library(tidyverse) # Viele Funktionen, u.a. für data.frames und ggplot2 für schöne plots
library(magrittr) # der Pipe-operator %>%: z.B: x = a; y=f(x); z=g(y); wird zu z= a %>% f %>% g
library(conveniencefunctions)
library(libSBML)
load("workspace.rda")Model setup
Read the model
This is basically copied from a libSBML-example
filename <- "/home/denial/Promotion/Projects/methacetin_fitting/model/met13_pkpd_7.xml" # Attention: tilde expansion doesn't work
d <- readSBML(filename)
m <- SBMLDocument_getModel(d)
# errors = SBMLDocument_getNumErrors(d);
# SBMLDocument_printErrors(d);
level = SBase_getLevel (d);
version = SBase_getVersion(d);
cat("File: ",filename," (Level ",level,", version ",version,")\n");
cat(" model id: ", ifelse(Model_isSetId(m), Model_getId(m) ,"(empty)"),"\n");
cat( "functionDefinitions: ", Model_getNumFunctionDefinitions(m) ,"\n" );
cat( " unitDefinitions: ", Model_getNumUnitDefinitions (m) ,"\n" );
cat( " compartmentTypes: ", Model_getNumCompartmentTypes (m) ,"\n" );
cat( " specieTypes: ", Model_getNumSpeciesTypes (m) ,"\n" );
cat( " compartments: ", Model_getNumCompartments (m) ,"\n" );
cat( " species: ", Model_getNumSpecies (m) ,"\n" );
cat( " parameters: ", Model_getNumParameters (m) ,"\n" );
cat( " initialAssignments: ", Model_getNumInitialAssignments (m) ,"\n" );
cat( " rules: ", Model_getNumRules (m) ,"\n" );
cat( " constraints: ", Model_getNumConstraints (m) ,"\n" );
cat( " reactions: ", Model_getNumReactions (m) ,"\n" );
cat( " events: ", Model_getNumEvents (m) ,"\n" );
cat( "\n" );Einschub: der Pipe-Operator %>%
# Mit dem Pipe-Operator kann man Funktionen verketten
# Standardmäßig wird das vorherige Ergebnis als erstes Argument von der nächsten Funktion eingesetzt. wenn man das nicht will, kann man es als . woanders hinsetzen
f <- function(x) x^2;
g <- function(x,y) x-y;
2 %>% f # f(2)
2 %>% f %>% g(3) # g(f(2),3) = 4-3 = 1
2 %>% f %>% g(3,.) # g(3,f(2)) = 3-4 = -1
# Man kann auch Funktionen definieren, die mit . losgehen, was dann das Argument für die Funktion ist.
h <- . %>% sqrt %>% add(5)
h(4)
# Das ist besonders nützlich in lapply, sapply und so weiter, wo man über eine liste(oder einen vektor) immer wieder die gleiche funktion laufen lässt (syntactic sugar für for-loops)
sapply(1:4, function(i) i^2 +5)
sapply(1:4 , . %>% raise_to_power(2) %>% add(5) ) # dassselbeAssignment Rules
- Get the assignment rules as string
- Apply the assignment rules onto themselves with str_replace. This is for assigments that convert parameters into other parameters
# get rules
nrules <- Model_getNumRules(m)
lrules <- Model_getListOfRules(m)
rules <- structure(sapply(0:(nrules-1), . %>% ListOfRules_get(lrules,.) %>% Rule_getFormula),
names = sapply(0:(nrules-1), . %>% ListOfRules_get(lrules,.) %>% Rule_getId))
rulenames <- names(rules)
# "Cure" rules: Since I do a parameter trafo for the units, I don't want to have any unit conversions via some rules
# A "bad" rule would be eg "QC = CO*3600/100", since I take care of the Units later. Therefore, the rule should be only "QC = CO"
rules <- rules %>% str_replace_all(c("1000" = "1", "3600" = "1", "\\b60\\b" = 1)) %>% set_names(rulenames)
# Apply the rules onto themselves to insert parameter transformations
# Final goal is to have a named vector where
# names are the "inner" parameters that are used within the model
# values are functions of "outer" parameters that are fed into the model
# apply rules 1st time
rules <- paste0("(", rules, ")") %>% set_names(paste0("\\b", rulenames, "\\b"))
rules <- str_replace_all(rules, rules) %>% set_names(rulenames)
# apply rules 2nd time
rules <- paste0("(", rules, ")") %>% set_names(paste0("\\b", rulenames, "\\b"))
rules <- str_replace_all(rules, rules) %>% set_names(rulenames)
# check if any of the rules are functions of other rules
indices <- rules %>% sapply(. %>% str_detect(paste0("\\b",rulenames, "\\b")) %>% any)
rules %>% extract(indices) %>% sapply(. %>% str_detect(paste0("\\b",rulenames, "\\b")) %>% extract(rulenames,.)) # check works, none of the
# print(rules)
# getSymbols(rules) # These are the "outer" parametersReactions
# get reactions
nreactions <- Model_getNumReactions(m)
lreactions <- Model_getListOfReactions(m)
reactions <- lapply(0:(nreactions-1), function(i) {
reaction <- ListOfReactions_get(lreactions,i)
nreactants <- reaction %>% Reaction_getNumReactants()
if (nreactants > 0) {
lreactants <- reaction %>% Reaction_getListOfReactants()
myreactant_stoichiometries <- lapply(0:(nreactants-1), . %>% ListOfSpeciesReferences_get(lreactants,.) %>% SpeciesReference_getStoichiometry)
myreactant_IDs <- lapply(0:(nreactants-1), . %>% ListOfSpeciesReferences_get(lreactants,.) %>% SimpleSpeciesReference_getSpecies)
from = paste0(paste0("(", myreactant_stoichiometries, "*", myreactant_IDs, ")"), collapse = "+")
} else {
from = ""
}
nproducts <- reaction %>% Reaction_getNumProducts()
if (nproducts > 0) {
lproducts <- reaction %>% Reaction_getListOfProducts()
myproduct_stoichiometries <- lapply(0:(nproducts-1), . %>% ListOfSpeciesReferences_get(lproducts,.) %>% SpeciesReference_getStoichiometry)
myproduct_IDs <- lapply(0:(nproducts-1), . %>% ListOfSpeciesReferences_get(lproducts,.) %>% SimpleSpeciesReference_getSpecies)
to = paste0(paste0("(", myproduct_stoichiometries, "*", myproduct_IDs, ")"), collapse = "+")
} else {
to = ""
}
# Apply rules to rate expressions
myrules <- rules
mynames <- names(myrules)
absorption_indices <- str_detect(myrules, "Absorption")
absorption_rules <- myrules[absorption_indices] %>% structure(names = mynames[absorption_indices])
myrules <- structure(myrules[!str_detect(myrules, "Absorption")], names = paste0("\\b", mynames[!str_detect(myrules, "Absorption")], "\\b"))
rate <- reaction %>% Reaction_getKineticLaw() %>% KineticLaw_getFormula() %>% str_replace_all(myrules)
description <- reaction %>% Reaction_getId()
# Incorporate the absorption:
# For this I add another reaction which absorbs the oral dose, e.g. D_apap_sul
if(description %>% str_detect("Absorption")) {
from[2] <- names(absorption_rules)[absorption_rules %>% str_detect(paste0(description, "\\)*$"))]
to[2] <- ""
rate[2] <- rate
description[2] <- description
}
# print(i)
return(data.frame(from = from, to = to, rate = rate, description = description, stringsAsFactors = F))
}) %>% do.call(rbind,.) # make one big data.frame out of it
# Build the dMod-object "eqnlist", which stores the reactions and stoichiometries
el <- eqnlist()
for(i in 1:nrow(reactions)) el <- addReaction(el, reactions$from[i], reactions$to[i], reactions$rate[i], reactions$description[i])
# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f <- el %>% as.eqnvec()
# This is the ODE when every rule is applied
print(f %>% head)Parameters
# get the parameters from the definition
nsbml_parameters <- m %>% Model_getNumParameters()
lsbml_parameters <- m %>% Model_getListOfParameters()
sbml_parameters <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getValue),
names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId))
sbml_parameter_descriptions <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getName),
names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId)) %>% {.[order(names(.))]}
# unit conversion
sbml_parameter_units <- structure(sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getUnits), names = sapply(0:(nsbml_parameters-1), . %>% ListOfParameters_get(lsbml_parameters,.) %>% Parameter_getId))
# The factors to bring each parameter to the units seconds, grams, litres
multiplication_factors <- sapply(0:(nsbml_parameters-1), . %>%
ListOfParameters_get(lsbml_parameters,.) %>%
Parameter_getDerivedUnitDefinition %>%
{myunitdef <- .
nunits <- UnitDefinition_getNumUnits(myunitdef)
lunits <- UnitDefinition_getListOfUnits(myunitdef)
lapply(0:(nunits-1), . %>%
ListOfUnits_get(lunits,.) %>%
{(Unit_getMultiplier(.) * 10^(Unit_getScale(.)))^Unit_getExponent(.)}) %>% Reduce("*",.)
} %>%
{.}
)
sbml_parameters <- sbml_parameters * multiplication_factors
# Initial assignments
niass <- m %>% Model_getNumInitialAssignments()
liass <- m %>% Model_getListOfInitialAssignments()
iass <- structure(sapply(0:(niass-1), . %>% ListOfInitialAssignments_get(liass,.) %>% print %>% str_split("\n") %>% sapply( . %>% {.[str_detect(.,"<ci>")]}) %>% str_replace_all("^.*<ci> ", "") %>% str_replace_all(" </ci>.*$", "")), names = sapply(0:(niass-1), . %>% ListOfInitialAssignments_get(liass,.) %>% InitialAssignment_getId))
nspecies <- Model_getNumSpecies(m)
lspecies <- Model_getListOfSpecies(m)
species <- lapply(0:(nspecies-1), . %>% {ListOfSpecies_get(lspecies,.)} %>% Species_getId()) %>% do.call(c,.)
inits <- structure(sapply(0:(nspecies-1), . %>% ListOfSpecies_get(lspecies,.) %>% Species_getInitialAmount), names = species)
inits[names(iass)] <- sbml_parameters[iass]
# All parameters combined, more than actually needed, because many of them have been replaced due to Rules
all_pars <- c(sbml_parameters, inits) # all possible parameters
# To check if the unit conversion worked
data.frame(par = names(sbml_parameters), value = sbml_parameters, original_unit= sbml_parameter_units, multiplication_factor = multiplication_factors, original_value = sbml_parameters/multiplication_factors, row.names = NULL) %>%
head()
#%>%
# filter(par %in% c("MET2APAP_HLM_CL", "fumic_metc13", "MPPGL", "BW", "FVli")) #%>% summarise(prod(multiplication_factor))
# filter(par %in% c("CO", "QC"))
# rules[str_detect(rules, "MET2APAP")]
# rules[str_detect(rules, "CO")]C-Code
the odemodel() command takes a while because sensitivity equations are calculated for derivatives and then the whole system is compiled into c-code.
# myodemodel <- odemodel(f, modelname = "methacetin") # this compiles the ode into a c-file. you can comment out this and the next line after it has been run once to save time.
# save(myodemodel, file = "methacetin.rda")
load("methacetin.rda")Prediction function
Make a prediction function from the odemodel. x will be a function x(times, pars)
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]Example plot
This plot is supposed to be the first plot from chunk 3 of the html-file that you sent me via email.
pars["Ave_metc13"] <- 0
pars["D_metc13"] <- 0.1 #100mg
mytimes = seq(0,8*3600, length.out = 80) #8 hours
# pred <- (g*x)(0:50, mypars, deriv = F, conditions = rownames(attr(data, "cond")))
prediction <- (x)(times = mytimes, pars = pars, deriv = F, conditions = "c1")
# Compute the volumes in litres for each state
volumes <- lapply(0:(nspecies-1), . %>% {ListOfSpecies_get(lspecies,.)} %>% Species_getCompartment()) %>% do.call(c,.) %>% set_names(species)
myrules <- rules %>% set_names(paste0("\\b", rulenames, "\\b"))
volumes <- str_replace_all(volumes, myrules) %>% set_names(species)
vol_fun <- funC0(volumes)
vol <- do.call(vol_fun, as.list(pars)) %>% t %>% {data.frame(volume=., name = rownames(.), stringsAsFactors = F)}
# Plot
plot_fun <- function(pred) pred %>%
wide2long() %>%
full_join(vol,by= "name") %>%
mutate(value = value * 1000, time = time/3600, concentration = value/volume) %>% # scale g to mg, s to h
separate(col = name, into = c("compartment", "substance"), sep = "_", extra = "merge") %>%
filter(compartment %in% c("Ali", "Ave", "Alu")) %>% # plot only liver, venuous and lung (as in the html file)
ggplot(aes(x= time, y = concentration)) +
geom_line(aes(color = compartment)) +
facet_wrap(~substance, scales = "free") +
theme_dMod() + scale_color_dMod()
plot_fun(prediction)
prediction_long <- (x)(times = mytimes*20,pars = pars, deriv = F, conditions = "c1")
plot_fun(prediction_long)
# Second plot from chunk 3
pars <- all_pars[getParameters(x)]
pars["Ave_metc13"] <- 0.1
pars["D_metc13"] <- 0 #100mg
# pars["Kplu_metc13"] <- 1
# pars["FVlu"] <- 7.6*10^(-6) * 10
prediction_second_plot <- (x)(times = mytimes,pars = pars, deriv = F, conditions = "c1")
plot_fun(prediction_second_plot)Modelling
Reduce the model complexity by inserting the fixed parameter values
free_parameters <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km" # Km value
)
fixed_parameters <- pars[!(names(pars)%in%c(free_parameters,names(f)[1]))] %>% names
myodemodel <- odemodel(f, modelname = "methacetin_reduced", fixed = fixed_parameters)
x <- Xs(myodemodel)Observation function
observables <- c(apap = "Ave_apap/(BW*FVve)",
apap_glu = "Ave_apap_glu/(BW*FVve)",
apap_sul = "Ave_apap_sul/(BW*FVve)",
apap_cys = "Ave_apap_cys/(BW*FVve)"
)
g <- Y(observables, x, parameters = free_parameters)(g*x)(mytimes, pars) %>% set_names("condition1") %>% plotPrediction(name %in% names(observables))Data
Load the data and transform it
myfiles <- list.files("~/Promotion/Projects/methacetin_fitting/data/", full.names = T)
raw_data <- myfiles %>% lapply(. %>% read.table(header = T, sep = "\t", stringsAsFactors = F))
data <-
raw_data %>%
lapply(. %>%
select(-contains("_mol")) %>%
gather("name_std", "std", ends_with("_sd")) %>%
mutate(name_std = str_replace(name_std, "_sd","")) %>%
gather("name_sigma", "sigma", ends_with("_se")) %>%
mutate(name_sigma = str_replace(name_sigma, "_se","")) %>%
{gather(.,"name", "value", one_of(.$name_std))} %>%
filter(name == name_std, name == name_sigma) %>%
# select(name,time,value,sigma) %>%
{.}) %>%
do.call(dMod::combine,.) %>%
mutate(D_apap = "D_apap", Ave_apap = "Ave_apap" ) %>%
{.$D_apap[.$study=="Chan1997"] <- 1400 / 1000
.$Ave_apap[.$study=="Chan1997"] <- 0
.$D_apap[.$study=="Chiew2010"] <- 5600 / 1000
.$Ave_apap[.$study=="Chiew2010"] <- 0
.$D_apap[.$study=="Critchley2005"] <- 1400 /1000
.$Ave_apap[.$study=="Critchley2005"] <- 0
.$D_apap[.$study=="Rawlins1977"] <- .$dose[.$study=="Rawlins1977"] * (.$route[.$study=="Rawlins1977"]%in%"oral") / 1000
.$Ave_apap[.$study=="Rawlins1977"] <- .$dose[.$study=="Rawlins1977"] * (.$route[.$study=="Rawlins1977"]%in%"iv") / 1000
.
} %>%
mutate(time = time * 3600, value = value/1000, sigma = sigma/1000) %>%
select(-group, -health_status, - name_std, - name_sigma, -std, -ethnicity, -route, -dose, -substance) %>%
# filter(!is.na(sigma)) %>%
# as.datalist() %>%
{.}
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
plot(mydatalist)Parameter transformations to define the conditions
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters] <- paste0("exp(log", free_parameters, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p <- NULL
for(i in 1:length(p_list)) { p <<- p + p_list[[i]]}
pouter <- log(pars[free_parameters]) %>% set_names(paste0("log",names(.)))
mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
plotCombined(mypred, mydatalist, name%in% names(observables))Fitting
obj <- normL2(mydatalist, (g*x*p))
# myfit <- mstrust(objfun = obj, center = pouter, studyname = "methacetin", cores = 3)
# save(myfit, file = "fit.rda")load("fit.rda")
fitted_pars <- myfit %>% as.parframe() %>% as.parvec()
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somypred <- (g*x*p)(seq(0, 48*3600, length.out = 100), fitted_pars)
liver <- names(f)[str_detect(names(f), "li")]
medication <- c(names(f)[str_detect(names(f), "D_apap$")], "Ave_apap")
myfit %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
myfit %>% as.parframe() %>% plotValues()
plotCombined(mypred, mydatalist, name %in% c(names(observables), liver, medication))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 8 2407.821 TRUE 60 -5.125968 -2.129845
## 2 35 2407.821 TRUE 71 -5.125973 -2.129851
## 3 21 2407.821 TRUE 39 -5.125965 -2.129842
## 4 29 2407.821 TRUE 67 -5.125962 -2.129839
## 5 33 2407.821 TRUE 76 -5.125966 -2.129843
## 6 22 2407.821 TRUE 75 -5.125969 -2.129847
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 -6.118641 -10.49793 -6.301154
## 2 -6.118641 -10.49793 -6.301161
## 3 -6.118641 -10.49794 -6.301223
## 4 -6.118641 -10.49794 -6.301221
## 5 -6.118641 -10.49794 -6.301214
## 6 -6.118641 -10.49794 -6.301198
## logAPAPCYS_HLM_CL: 2.75936204923354e-05
## logAPAPCYS_Km: 0.00183418625326802
## logAPAPGLU_HLM_CL: 0.00594046463001759
## logAPAPGLU_Km: 0.118855661904664
## logAPAPSUL_HLM_CL: 0.00220144570190944Free some additional parameters
Free other parameters try 1 - not good
Here I freed some additional parameters, but since I freed “Ka_apap_sul” instead of “Ka_apap”, the results don’t make much additional sense.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters1 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul", "F_apap_sul",
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu",
"Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
"Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters1 <- pars[!(names(pars)%in%c(free_parameters1,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters1] <- paste0("exp(log", free_parameters1, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p1 <- NULL
for(i in 1:length(p_list)) { p1 <<- p1 + p_list[[i]]}
pouter <- log(pars[free_parameters1]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj1 <- normL2(mydatalist, (g*x*p1))
# job1 <- runbg({myfit <- mstrust(objfun = obj1, center = pouter, studyname = "methacetin", cores = 20, fits = 100); myfit}, machine = "knecht3", filename = "job1")# myfit1 <- job1$get()$knecht3
# save(myfit1, file = "myfit1.rda")
# job1$purge()
load("myfit1.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit1 %>% as.parframe())+scale_y_log10()
mypred1 <- (g*x*p1)(mytimes*4, myfit1 %>% as.parframe() %>% as.parvec)
plotCombined(mypred1, mydatalist, name %in% names(observables))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 71 2749.704 FALSE 100 -3.121425 0.9276677
## 2 17 2871.109 FALSE 100 -4.825791 -0.8966849
## 3 4 2879.499 TRUE 82 -4.836514 -2.0500243
## 4 83 3006.360 TRUE 38 -4.810375 -1.6071039
## 5 11 3025.134 TRUE 47 -5.068367 -2.1571283
## 6 60 3063.915 FALSE 100 -4.778175 -1.9189665
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -5.966762 -8.877054 -8.056205 -9.043602
## 2 -7.312292 -9.776533 -6.317704 -7.411144
## 3 -6.424722 -9.745092 -7.247758 -7.449345
## 4 -6.432805 -10.767499 -6.451981 -6.734894
## 5 -6.368207 -9.236633 -7.298571 -8.975554
## 6 -6.479809 -10.151009 -5.692777 -5.931667
## logF_apap_sul logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys
## 1 0.53130523 -0.2895295 -0.9819923 0.32392147 3.4721895
## 2 0.07956807 -0.3510074 0.2160174 0.99467641 0.7860616
## 3 -0.52122381 -0.3932631 -1.0615523 -0.09235355 0.9514369
## 4 0.16252315 -0.3785620 0.1383524 0.27175050 -0.8202826
## 5 1.21332440 -0.2100422 -1.7021606 -0.01216570 1.9565013
## 6 -2.38273975 -0.4474130 -1.6230284 0.04538435 0.4837584
## logKpki_apap_cys logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu
## 1 -0.9981306 -0.7078988 -1.1713447 -0.4586217
## 2 1.4876644 0.5391938 0.5938409 0.6894229
## 3 1.3055012 0.4413984 0.7986194 0.3810163
## 4 -0.1222928 -0.2784549 -0.6756876 0.1384249
## 5 1.2069110 0.6198618 0.5661661 -1.1336736
## 6 -0.1000206 -0.1492559 0.6852967 -0.2904269
## logKpli_apap_glu logKpre_apap_sul logKpre_co2c13 logKpli_co2c13
## 1 1.2820601 0.8610661 -1.8050599 -1.0958784
## 2 -1.2743622 -2.0279402 -1.6334736 0.2964431
## 3 -0.7746656 -1.5237195 0.3130873 0.3634344
## 4 0.4949366 -0.5271704 -0.8464643 -1.1150170
## 5 -1.4231745 -1.4809556 0.9901271 -0.4769397
## 6 1.7057846 -1.3223539 -1.2621168 -0.9185781
## logKpre_metc13 logKpli_metc13
## 1 -2.4242875 -0.4355059
## 2 -1.3074538 -0.5878608
## 3 0.0672823 1.2834798
## 4 -2.5253310 1.1213766
## 5 -1.0586655 1.3624494
## 6 -1.3302118 1.4817382
## logAPAPCYS_HLM_CL: 0.000139554700132717
## logAPAPCYS_Km: 0.000317128135451512
## logAPAPGLU_HLM_CL: 0.0440942903336689
## logAPAPGLU_Km: 2.52860494207782
## logAPAPSUL_HLM_CL: 0.00256252477737589
## logF_apap_sul: 1.70115125864559
## logKa_apap_sul: 0.000118144550067709
## logKpki_apap: 0.374564129323335
## logKpki_apap_cys: 0.368567801198445
## logKpki_apap_glu: 0.63215432668864
## logKpli_apap: 1.38253873675798
## logKpli_apap_cys: 0.492678343277563
## logKpli_apap_glu: 3.60405672997012
## logKpli_co2c13: 0.334245860035197
## logKpli_metc13: 0.646937302760331
## logKpre_apap: 0.748615739138366
## logKpre_apap_cys: 32.2071824118807
## logKpre_apap_glu: 0.309949884477503
## logKpre_apap_sul: 2.36568138653131
## logKpre_co2c13: 0.164464609213266
## logKpre_metc13: 0.0885411804880305Free other parameters 2 - not good
Here I freed some less parameters than in try 1 but it suffers from the same problem
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters2 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul", "F_apap_sul" #, # total dumm, Ka und F sind so redundant wie es nur geht.
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu",
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters2 <- pars[!(names(pars)%in%c(free_parameters2,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters2] <- paste0("exp(log", free_parameters2, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p2 <- NULL
for(i in 1:length(p_list)) { p2 <<- p2 + p_list[[i]]}
pouter <- log(pars[free_parameters2]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj2 <- normL2(mydatalist, (g*x*p2))
# job2 <- runbg({myfit <- mstrust(objfun = obj2, center = pouter, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht4", filename = "job2")
# # myfit2 <- job2$get()$knecht4
# save(myfit2, file = "myfit2.rda")
# job2$purge()
load("myfit2.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit2 %>% as.parframe())+scale_y_log10()
mypred2 <- (g*x*p2)(mytimes*4, myfit2 %>% as.parframe() %>% as.parvec)
plotCombined(mypred2, mydatalist, name %in% names(observables))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 24 3134.039 TRUE 25 -5.067594 -2.064909
## 2 56 3134.040 TRUE 30 -5.067589 -2.064904
## 3 20 3134.040 TRUE 26 -5.067595 -2.064908
## 4 89 3134.040 TRUE 16 -5.067592 -2.064905
## 5 37 3134.040 TRUE 86 -5.067589 -2.064901
## 6 32 3134.040 TRUE 20 -5.067590 -2.064902
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -6.117757 -10.48519 -6.241245 -7.346186
## 2 -6.117757 -10.48520 -6.241296 -7.021465
## 3 -6.117758 -10.48501 -6.240447 -7.216762
## 4 -6.117758 -10.48501 -6.240447 -6.306130
## 5 -6.117757 -10.48501 -6.240481 -7.986534
## 6 -6.117757 -10.48501 -6.240461 -7.899805
## logF_apap_sul
## 1 41.111322
## 2 43.933670
## 3 50.970546
## 4 43.726021
## 5 7.397659
## 6 34.039840
## logAPAPCYS_HLM_CL: 2.79472212506865e-05
## logAPAPCYS_Km: 0.00194743024961977
## logAPAPGLU_HLM_CL: 0.00629755070545838
## logAPAPGLU_Km: 0.126829778557011
## logAPAPSUL_HLM_CL: 0.00220339319347716
## logF_apap_sul: 715187883400265600
## logKa_apap_sul: 0.000645048206917656Free other parameters 3 - not good
Same problem as in tries 1 and 2
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters3 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul", "F_apap_sul" ,
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters3 <- pars[!(names(pars)%in%c(free_parameters3,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters3] <- paste0("exp(log", free_parameters3, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p3 <- NULL
for(i in 1:length(p_list)) { p3 <<- p3 + p_list[[i]]}
pouter <- log(pars[free_parameters3]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj3 <- normL2(mydatalist, (g*x*p3))
# job3 <- runbg({myfit <- mstrust(objfun = obj3, center = pouter, studyname = "methacetin", cores = 15, fits = 100); myfit}, machine = "knecht1", filename = "job3")# myfit3 <- job3$get()$knecht1
# save(myfit3, file = "myfit3.rda")
# job3$purge()
load("myfit3.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit3 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit3 %>% as.parframe())+scale_y_log10()
mypred3 <- (g*x*p3)(mytimes*4, myfit3 %>% as.parframe() %>% as.parvec)
plotCombined(mypred3, mydatalist, name %in% names(observables))
## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 59 2562.482 TRUE 84 -5.201098 -2.097660
## 2 60 2689.397 TRUE 48 -5.383774 -2.043510
## 3 55 2856.425 FALSE 100 -5.109883 -1.468816
## 4 48 3019.107 FALSE 100 -6.865873 -5.295360
## 5 1 3039.058 FALSE 100 -4.638935 -1.383292
## 6 89 3066.281 TRUE 76 -6.464966 -4.620926
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -6.555927 -10.709675 -6.33769669 -6.042436
## 2 -6.862922 -10.910624 -6.20766738 -8.289160
## 3 -6.656348 -9.872593 -4.75970985 -6.579562
## 4 -6.050377 -9.618928 -6.86802281 -6.542891
## 5 -6.709802 -7.521883 -0.09424782 -7.980098
## 6 -5.857182 -10.682944 -8.86895808 -7.943206
## logF_apap_sul logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys
## 1 2.808457812 -0.3669234 0.7613894 0.1672635 -0.6048195
## 2 -1.512335536 -0.1336978 -0.2052586 0.5449920 -1.5425344
## 3 0.009716211 -0.2655233 -1.6898116 0.6768234 -0.2769801
## 4 -0.661252544 -0.1334166 0.2229806 -0.6106942 -0.4045914
## 5 2.748915536 -0.4468966 -0.6003575 0.3837791 -0.3193860
## 6 1.679453017 -0.1192626 -0.5590351 -0.6042076 -2.9191368
## logKpki_apap_cys logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu
## 1 -0.1241948 -0.04614371 -0.8681152 0.36911867
## 2 -0.2511150 1.52077842 -0.8083817 0.16712529
## 3 0.7221217 -0.09379088 -1.0941684 0.26637337
## 4 2.2752818 0.30050894 -1.5193795 1.64286021
## 5 -2.5272957 1.04456393 0.3086699 0.04221941
## 6 0.4521719 -0.72613317 0.1994445 0.85658306
## logKpli_apap_glu logKpre_apap_sul
## 1 -0.003140001 -1.0894804
## 2 -0.342208285 -2.2243255
## 3 -0.207101279 -0.1223363
## 4 -0.702140533 -2.0416762
## 5 0.408290615 -1.1143709
## 6 -0.820070493 -0.7383285
## logAPAPCYS_HLM_CL: 2.23278620657747e-05
## logAPAPCYS_Km: 0.00176837065068179
## logAPAPGLU_HLM_CL: 0.00551051224903083
## logAPAPGLU_Km: 0.122743302983146
## logAPAPSUL_HLM_CL: 0.0014216640114061
## logF_apap_sul: 16.5843223443174
## logKa_apap_sul: 0.00237576415394654
## logKpki_apap: 2.14124916301601
## logKpki_apap_cys: 0.883207793863882
## logKpki_apap_glu: 1.44645924035911
## logKpli_apap: 1.18206573439193
## logKpli_apap_cys: 0.954904720644904
## logKpli_apap_glu: 0.996864923604179
## logKpre_apap: 0.692862701926673
## logKpre_apap_cys: 0.546172985593896
## logKpre_apap_glu: 0.419741938691713
## logKpre_apap_sul: 0.336391236317653Free other parameters 4 - not good
Same problem as in tries 1,2,3. The only difference to try 2 is that I removed the structural non-identifiability of Ka_apap_sul and F_apap_sul.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters4 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap_sul"#, "F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters4 <- pars[!(names(pars)%in%c(free_parameters4,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters4] <- paste0("exp(log", free_parameters4, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p4 <- NULL
for(i in 1:length(p_list)) { p4 <<- p4 + p_list[[i]]}
pouter <- log(pars[free_parameters4]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj4 <- normL2(mydatalist, (g*x*p4))
# job4 <- runbg({myfit <- mstrust(objfun = obj4, center = pouter, studyname = "methacetin", cores = 15, fits = 100); myfit}, machine = "knecht1", filename = "job4")
# save(job4, file = "job4.rda")# myfit4 <- job4$get()$knecht1
# save(myfit4, file = "myfit4.rda")
# job4$purge()
load("myfit4.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit4 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit4 %>% as.parframe())+scale_y_log10()
mypred4 <- (g*x*p4)(mytimes*4, myfit4 %>% as.parframe() %>% as.parvec, deriv = F)
myplot4 <- plotCombined(mypred4, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot4)
# myplot4## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 65 3134.039 TRUE 31 -5.067596 -2.064909
## 2 47 3134.040 TRUE 74 -5.067589 -2.064901
## 3 51 3134.040 TRUE 82 -5.067591 -2.064903
## 4 49 3134.040 TRUE 58 -5.067589 -2.064902
## 5 17 3134.040 TRUE 18 -5.067592 -2.064904
## 6 91 3134.040 TRUE 25 -5.067590 -2.064903
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap_sul
## 1 -6.117758 -10.48501 -6.240449 69.0940136
## 2 -6.117757 -10.48502 -6.240488 -1.9728623
## 3 -6.117757 -10.48501 -6.240450 -1.9723416
## 4 -6.117757 -10.48501 -6.240446 -0.7434353
## 5 -6.117757 -10.48501 -6.240455 21.5604671
## 6 -6.117758 -10.48499 -6.240395 39.2232311
## logAPAPCYS_HLM_CL: 2.79524088116887e-05
## logAPAPCYS_Km: 0.00194898079879723
## logAPAPGLU_HLM_CL: 0.0062975431281648
## logAPAPGLU_Km: 0.126829823077826
## logAPAPSUL_HLM_CL: 0.00220339066957284
## logKa_apap_sul: 1.01659705294599e+30Further hypotheses
Freeing the parameters in the upper four fits didn’t help much. On the other hand, lots of the fits didn’t converge. From here, I could do several things.
- Introduce scaling factors for each condition
- Fill the empty sigmas with fitted sigmas
- Run a profile likelihood analysis for the fits 1-4
- Re-run the fits 2 and 3 without the structural non-identifiability of Ka_apap_sul and F_apap_sul
- If I re-run the fits, I can take the fitted values of the original free parameters as center for the sampling
- In data, drop the column n before converting it to a datalist
- In the upper fits, I freed the parameters Ka_apap_sul, which of course doesn’t make sense, since no apap_sul is given. I have to free the parameters Ka_apap
Free other parameters 5 - not good
Here I freed Ka_apap which makes sense conceptually (increase absorption rate). But it doesn’t improve the fits much.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters5 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, "F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters5 <- pars[!(names(pars)%in%c(free_parameters5,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters5] <- paste0("exp(log", free_parameters5, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p5 <- NULL
for(i in 1:length(p_list)) { p5 <<- p5 + p_list[[i]]}
pouter <- log(pars[free_parameters5]) %>% set_names(paste0("log",names(.)))
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj5 <- normL2(mydatalist, (g*x*p5))
# job5 <- runbg({myfit <- mstrust(objfun = obj5, center = pouter, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht1", filename = "job5")
# save(job5, file = "job5.rda")# myfit5 <- job5$get()$knecht1
# save(myfit5, file = "myfit5.rda")
# job5$purge()
load("myfit5.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit5 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit5 %>% as.parframe())+scale_y_log10()
mypred5 <- (g*x*p5)(mytimes*5, myfit5 %>% as.parframe() %>% as.parvec, deriv = F)
myplot5 <- plotCombined(mypred5, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot5)
# myplot5## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 33 3092.015 FALSE 100 -5.048048 -2.013030
## 2 12 3092.025 TRUE 51 -5.047014 -2.013058
## 3 99 3092.031 FALSE 100 -5.048425 -2.014686
## 4 60 3092.117 TRUE 96 -5.047379 -2.011610
## 5 88 3092.146 TRUE 49 -5.047300 -2.011395
## 6 26 3092.203 TRUE 40 -5.048943 -2.016184
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -6.142913 -10.48392 -6.226134 -6.853661
## 2 -6.141722 -10.48302 -6.222827 -6.879080
## 3 -6.141694 -10.48360 -6.225265 -6.880146
## 4 -6.143583 -10.48343 -6.224064 -6.838579
## 5 -6.143710 -10.48332 -6.223574 -6.835689
## 6 -6.140809 -10.48370 -6.225923 -6.898840
## logAPAPCYS_HLM_CL: 2.79827234659552e-05
## logAPAPCYS_Km: 0.00197708080993727
## logAPAPGLU_HLM_CL: 0.00642185876478594
## logAPAPGLU_Km: 0.133583353550372
## logAPAPSUL_HLM_CL: 0.0021486564154929
## logKa_apap: 0.00105558386269278Free other parameters 6 - not enough fits converged
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters6 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap", #"F_apap_sul" ,
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters6 <- pars[!(names(pars)%in%c(free_parameters6,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters6] <- paste0("exp(log", free_parameters6, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p6 <- NULL
for(i in 1:length(p_list)) { p6 <<- p6 + p_list[[i]]}
pouter <- log(pars[free_parameters6]) %>% set_names(paste0("log",names(.)))
best_fit <- myfit %>% as.parframe() %>% as.parvec()
pouter[names(best_fit)] <- best_fit
# mypred <- (g*x*p)(seq(0, 48*3600, length.out = 200), pouter)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj6 <- normL2(mydatalist, (g*x*p6))
# job6 <- runbg({myfit <- mstrust(objfun = obj6, center = pouter, studyname = "methacetin", cores = 10, fits = 100); myfit}, machine = "knecht3", filename = "job6")
# save(job6, file = "job6.rda")
# job6$check()# myfit6 <- job6$get()$knecht3
# save(myfit6, file = "myfit6.rda")
# job6$purge()
load("myfit6.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit6 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit6 %>% as.parframe())+scale_y_log10()
mypred6 <- (g*x*p6)(mytimes*4, myfit6 %>% as.parframe() %>% {.[2,]} %>% as.parvec)
myplot6 <- plotCombined(mypred6, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot6)
# myplot6## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 16 2654.322 TRUE 88 -6.288191 -4.0064053
## 2 60 3025.841 FALSE 100 5.036281 8.1988600
## 3 39 3059.748 FALSE 100 -3.427671 0.2086894
## 4 19 3210.587 TRUE 100 -4.251426 -1.6417543
## 5 47 3233.777 TRUE 48 -5.809100 -2.7455612
## 6 48 3302.297 FALSE 100 -4.770855 -0.8605885
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -6.242811 -9.413691 -6.941120 -5.245524
## 2 -5.281948 -10.686351 -6.408999 -7.197151
## 3 -6.146221 -9.736307 -7.455677 -6.582344
## 4 -4.910387 -10.540902 -7.410143 -7.443229
## 5 -6.315936 -10.030452 -1.701016 -7.739147
## 6 -7.215143 -9.338544 -1.877134 -7.204190
## logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys logKpki_apap_cys
## 1 -0.04690208 0.1358953 -0.3757922 2.19032710 -0.4659840
## 2 -0.27191200 1.4639032 -0.5538510 -0.52221741 -0.3912030
## 3 -0.42939025 0.3304541 0.2229900 0.05838160 1.8274448
## 4 -0.18759530 -1.7257917 -0.7659503 -0.40658360 0.1236754
## 5 -0.08418012 0.4028411 0.2491479 -3.17511619 -2.8534696
## 6 -0.34371554 -1.2734140 1.0683678 0.03096511 -1.5696406
## logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu logKpli_apap_glu
## 1 -0.298471419 -1.2802038 0.54299839 0.09378855
## 2 0.276192230 -1.6506727 -0.45872701 -0.11353331
## 3 0.279366410 -1.3550508 -0.31509066 -0.10951404
## 4 0.009838299 -0.8169902 -0.19239433 -0.32722341
## 5 -0.007563093 -2.0011674 0.04291508 -0.71993332
## 6 -0.671261233 0.5743076 -1.17050260 1.01329858
## logKpre_apap_sul
## 1 -1.50864561
## 2 0.48599584
## 3 -0.05666556
## 4 0.78554977
## 5 -0.54997601
## 6 -1.66492170
## logAPAPCYS_HLM_CL: 8.15992019166072e-05
## logAPAPCYS_Km: 0.000967185735063076
## logAPAPGLU_HLM_CL: 0.00185811908010114
## logAPAPGLU_Km: 0.0181986967501468
## logAPAPSUL_HLM_CL: 0.00194438305778717
## logKa_apap: 0.0052710585929262
## logKpki_apap: 1.14556190846099
## logKpki_apap_cys: 0.627517333300021
## logKpki_apap_glu: 1.72115983885983
## logKpli_apap: 0.686745011371109
## logKpli_apap_cys: 0.741951487277269
## logKpli_apap_glu: 1.09832748267348
## logKpre_apap: 0.954180830717175
## logKpre_apap_cys: 8.93813625877344
## logKpre_apap_glu: 0.277980630583032
## logKpre_apap_sul: 0.22120937811377Free other parameters 6_1 - more iterations didn’t help, it’s still bad
In comparison to try 6, I also drop column n before converting the data to a datalist. I use the best fit from try 6 as center for the sampling
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters6 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap", #"F_apap_sul" ,
"Kpre_apap", "Kpki_apap", "Kpli_apap",
"Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
"Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
"Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters6 <- pars[!(names(pars)%in%c(free_parameters6,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
p_list <- lapply(1:nrow(conditions), function(i) {
trafo <- as.character(pars) %>% set_names(names(pars))
cond <- unlist(conditions[i,])[2:3]
trafo[names(cond)] <- cond
trafo[free_parameters6] <- paste0("exp(log", free_parameters6, ")")
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p6_1 <- NULL
for(i in 1:length(p_list)) { p6_1 <<- p6_1 + p_list[[i]]}
pouter6_1 <- log(pars[free_parameters6]) %>% set_names(paste0("log",names(.)))
best_fit <- myfit6 %>% as.parframe() %>% {.[2,]} %>% as.parvec() # der 2. sieht viel besser aus von der dynamik her
pouter6_1[names(best_fit)] <- best_fit
# mypred <- (g*x*p6_1)(seq(0, 48*3600, length.out = 200), pouter6_1, deriv = F)
# plotCombined(mypred, mydatalist, name%in% names(observables))
obj6_1 <- normL2(mydatalist, (g*x*p6_1))
# job6_1 <- runbg({myfit <- mstrust(objfun = obj6_1, center = pouter6_1, studyname = "methacetin", cores = 12, fits = 200, iterlim = 200, sd = 2); myfit}, machine = "knecht5", filename = "job6_1", input = global_env_without(c("mypred", "job", "myfit")))
# save(job6_1, file = "job6_1.rda")
# job6_1$check()# load("job6_1.rda")
# myfit6_1 <- job6_1$get()$knecht5
# save(myfit6_1, file = "myfit6_1.rda")
# job6_1$purge()
load("myfit6_1.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit6_1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit6_1 %>% as.parframe())+scale_y_log10()
mypred6_1 <- (g*x*p6_1)(mytimes*4, myfit6_1 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot6_1 <- plotCombined(mypred6_1, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot6_1)
# myplot6_1## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 126 2254.796 TRUE 168 -7.160945 -4.394436
## 2 55 2616.813 FALSE 200 -6.451695 -7.616642
## 3 17 2724.360 TRUE 107 5.557814 6.673239
## 4 40 3089.612 TRUE 86 4.825572 9.344937
## 5 49 3140.781 FALSE 200 -7.690501 -8.975971
## 6 152 3143.196 TRUE 57 -7.021284 -1.557577
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -6.718847 -8.311169 -1.954417 -6.524918
## 2 -2.815882 -9.898235 -11.899718 -4.243046
## 3 -4.036551 -9.848632 -6.435657 -7.323644
## 4 -5.141421 -9.827778 -10.177033 -6.097087
## 5 -4.286810 -7.639483 -2.995771 -6.741065
## 6 -7.080733 -10.818246 -6.381171 -5.752037
## logKpre_apap logKpki_apap logKpli_apap logKpre_apap_cys logKpki_apap_cys
## 1 0.011772621 1.1690414 0.4379879 0.9132539 -1.9017606
## 2 -0.217629833 1.1177705 -3.9170836 0.4071650 1.7266591
## 3 -0.270321088 -1.2532600 -2.0291797 0.6635614 -1.2993619
## 4 -0.296894373 2.2177932 -0.7600494 -1.0904734 2.5744728
## 5 0.005178141 1.4576398 -2.2549729 -0.7537454 -2.2248627
## 6 0.234766940 0.9271649 1.8868622 -1.9757345 0.2620245
## logKpli_apap_cys logKpre_apap_glu logKpki_apap_glu logKpli_apap_glu
## 1 -3.2571647 -1.58858870 0.4221563 -0.4352780
## 2 0.0398897 -1.46330331 1.4233708 -0.4196139
## 3 -2.7804837 0.05884825 -0.1986642 -1.6136809
## 4 1.5410468 -3.05246767 -2.0834014 -4.2756599
## 5 2.8185400 -1.92875393 0.2224316 5.5323154
## 6 -1.5929422 -3.93035938 -0.3190701 -3.6514846
## logKpre_apap_sul
## 1 -0.6014329
## 2 -0.4671365
## 3 0.2192896
## 4 1.7863395
## 5 -0.4661242
## 6 1.1538600
## logAPAPCYS_HLM_CL: 0.000245756553408302
## logAPAPCYS_Km: 0.141647012033554
## logAPAPGLU_HLM_CL: 0.000776320239293284
## logAPAPGLU_Km: 0.0123458440700416
## logAPAPSUL_HLM_CL: 0.00120793042044707
## logKa_apap: 0.0014664396901935
## logKpki_apap: 3.21890543991051
## logKpki_apap_cys: 0.1493055246381
## logKpki_apap_glu: 1.52524686894739
## logKpli_apap: 1.54958613737682
## logKpli_apap_cys: 0.038497394575441
## logKpli_apap_glu: 0.647084766249551
## logKpre_apap: 1.01184219056339
## logKpre_apap_cys: 2.4924194943794
## logKpre_apap_glu: 0.204213614409576
## logKpre_apap_sul: 0.548025794974965global_env_without <- function(reg) ls(.GlobalEnv)[!(ls(.GlobalEnv) %>% sapply(. %>% str_detect(reg) %>% any))]Scaling factors and Error model
Introduce scaling factors 7 - not enough fits converged, but in principle not bad
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters7 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters7 <- pars[!(names(pars)%in%c(free_parameters7,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables7 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters7 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters7 <- c(free_parameters7, scale_parameters7)
i <- 2
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters7] <- paste0("exp(log", free_parameters7, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters7, x = scale_parameters7, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters7)
scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any %>% not)] <- "1"
trafo <- c(trafo, scales)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p7 <- NULL
for(i in 1:length(p_list)) { p7 <<- p7 + p_list[[i]]}
g7 <- Y(observables7, x)#, parameters = c(free_parameters7, scale_parameters7))
obj7 <- normL2(mydatalist, (g7*x*p7))
pouter7 <- rep(0, length(getParameters(obj7))) %>% set_names(getParameters(obj7))
pouter7[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# job7 <- runbg({myfit <- mstrust(objfun = obj7, center = pouter7, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht5", filename = "job7", input = global_env_without(c("mypred", "job", "myfit")))
# save(job7, file = "job7.rda")
# job7$check()# myfit7 <- job7$get()$knecht5
# save(myfit7, file = "myfit7.rda")
# job7$purge()
load("myfit7.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit7 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit7 %>% as.parframe())+scale_y_log10()
mypred7 <- (g7*x*p7)(mytimes*4, myfit7 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot7 <- plotCombined(mypred7, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot7)
# myplot7## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 5 2490.012 TRUE 30 -5.749303 -3.08097475
## 2 64 2626.769 TRUE 60 -4.454476 -2.31516013
## 3 67 2855.601 TRUE 87 -5.293002 -2.39553148
## 4 35 3037.778 FALSE 100 -4.055656 -2.26447288
## 5 98 3347.208 FALSE 100 -5.763790 -4.18377196
## 6 14 3572.777 FALSE 100 -2.353277 -0.00949594
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -5.069446 -10.95827 -6.696778 -7.438818
## 2 -5.574571 -12.31058 -6.987258 -7.603783
## 3 -5.345860 -11.73765 -6.328117 -6.926927
## 4 -5.365076 -11.24999 -7.009670 -7.954771
## 5 -5.881820 -10.51708 -7.158013 -8.062787
## 6 -5.624627 -10.82433 -6.874139 -7.609710
## log_scale_apap_Chan1997_1.4_0 log_scale_apap_glu_Chan1997_1.4_0
## 1 0.7357347 -0.258980454
## 2 0.9075704 0.003338643
## 3 0.4662637 0.365011207
## 4 1.2621657 -0.338130480
## 5 1.0203440 -0.713985018
## 6 0.7625238 0.359423877
## log_scale_apap_sul_Chan1997_1.4_0 log_scale_apap_cys_Chan1997_1.4_0
## 1 -0.03917457 0.6911870
## 2 0.21839269 0.4441861
## 3 -0.26185349 1.4962657
## 4 0.48438243 0.9564866
## 5 0.35501653 0.2993361
## 6 -0.08656187 0.4475230
## log_scale_apap_Chiew2010_5.6_0 log_scale_apap_glu_Chiew2010_5.6_0
## 1 0.28907910 1.4721450
## 2 0.32685873 0.8062451
## 3 0.04167644 0.5089211
## 4 0.28628828 0.1863739
## 5 -0.25747805 -0.5369452
## 6 0.60731079 -1.4552898
## log_scale_apap_sul_Chiew2010_5.6_0 log_scale_apap_Critchley2005_1.4_0
## 1 -0.1620631 -0.3123042
## 2 0.8175662 0.1778291
## 3 -2.1840586 -0.8736844
## 4 0.8380141 0.2651291
## 5 0.2400412 -1.6486526
## 6 0.6459171 0.1341244
## log_scale_apap_glu_Critchley2005_1.4_0
## 1 0.6340540
## 2 -0.8673394
## 3 -0.5829681
## 4 0.1720526
## 5 -0.6356093
## 6 -1.1555520
## log_scale_apap_sul_Critchley2005_1.4_0
## 1 -0.32952950
## 2 -0.58790653
## 3 -0.75363120
## 4 -0.42899409
## 5 0.36572246
## 6 0.06592612
## log_scale_apap_cys_Critchley2005_1.4_0 log_scale_apap_Rawlins1977_0_1
## 1 0.4520727 0.12354755
## 2 1.9508658 0.19603601
## 3 1.1278751 0.01086081
## 4 0.8794284 0.34127665
## 5 -0.4761007 0.23294843
## 6 0.3483191 0.15601077
## log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 0.27999826 0.23205478
## 2 0.40194965 0.40054980
## 3 0.07005591 0.03746354
## 4 0.69124502 0.72483500
## 5 0.55520067 0.47151149
## 6 0.28847757 0.35586546
## log_scale_apap_Rawlins1977_0.5_0
## 1 -0.32097913
## 2 -0.03674230
## 3 -0.03173343
## 4 -0.78815698
## 5 0.17669530
## 6 0.45105346
## logAPAPCYS_HLM_CL: 1.74133804776165e-05
## logAPAPCYS_Km: 0.00123488444599134
## logAPAPGLU_HLM_CL: 0.0031850000915229
## logAPAPGLU_Km: 0.0459144794614278
## logAPAPSUL_HLM_CL: 0.00628590132812984
## logKa_apap: 0.000587979808605786
## log_scale_apap_Chan1997_1.4_0: 2.08701470180148
## log_scale_apap_Chiew2010_5.6_0: 1.3351973371447
## log_scale_apap_Critchley2005_1.4_0: 0.731758900597009
## log_scale_apap_cys_Chan1997_1.4_0: 1.99608345214162
## log_scale_apap_cys_Critchley2005_1.4_0: 1.57156625159931
## log_scale_apap_glu_Chan1997_1.4_0: 0.771838109420778
## log_scale_apap_glu_Chiew2010_5.6_0: 4.35857412292258
## log_scale_apap_glu_Critchley2005_1.4_0: 1.8852379216948
## log_scale_apap_Rawlins1977_0_1: 1.13150380172748
## log_scale_apap_Rawlins1977_0.5_0: 0.725438390393933
## log_scale_apap_Rawlins1977_1_0: 1.32312750633397
## log_scale_apap_Rawlins1977_2_0: 1.26118881349357
## log_scale_apap_sul_Chan1997_1.4_0: 0.96158283149328
## log_scale_apap_sul_Chiew2010_5.6_0: 0.850387509055339
## log_scale_apap_sul_Critchley2005_1.4_0: 0.719262064522545Error model
8 Introduce error model and scaling factors in the dynamic model - very weird results.
Very bad results… Mistake: Don’t remove the rows with sigma=NA! Stupid…
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters8 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters8 <- pars[!(names(pars)%in%c(free_parameters8,names(f)[1]))] %>% names
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables8 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters8 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters8 <- c(free_parameters8, scale_parameters8)
error_model8 <- c(apap = "srel_apap*apap^2 +s0_apap",
apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters8 <- setdiff(getSymbols(error_model8), names(error_model8)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters8] <- paste0("exp(log", free_parameters8, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters8, x = scale_parameters8, y = .)}
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters8, x = error_parameters8, y = .)}
errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p8 <- NULL
for(i in 1:length(p_list)) { p8 <<- p8 + p_list[[i]]}
g8 <- Y(observables8, x)#, parameters = c(free_parameters8, scale_parameters8))
err8 <- Y(error_model8, g8)
obj8 <- normL2(mydatalist, (g8*x*p8), errmodel = err8)
pouter8 <- rep(0, length(getParameters(obj8))) %>% set_names(getParameters(obj8))
pouter8[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# obj8(pouter8)
# job8 <- runbg({myfit <- mstrust(objfun = obj8, center = pouter8, studyname = "methacetin", cores = 12, fits = 100); myfit}, machine = "knecht5", filename = "job8", input = global_env_without(c("mypred", "job", "myfit")))
# save(job8, file = "job8.rda")
# job8$purge()# myfit8 <- job8$get()$knecht5
# save(myfit8, file = "myfit8.rda")
# job8$purge()
load("myfit8.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit8 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit8 %>% as.parframe())#+scale_y_log10()
mypred8 <- (g8*x*p8)(mytimes*4, myfit8 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot8 <- plotCombined(mypred8, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot8)
# myplot8## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 46 -486.5434 TRUE 68 0.8448040 3.0690138
## 2 95 -477.9974 TRUE 86 2.3221006 4.5447823
## 3 3 -458.5112 TRUE 85 -2.8788000 -0.5052364
## 4 91 -451.1060 TRUE 84 2.0359641 4.2579075
## 5 50 -437.6416 FALSE 100 2.1876617 4.4051112
## 6 44 -426.7063 TRUE 82 -0.9089178 1.3108703
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -10.788256 -13.56917 -3.126300 -6.295179
## 2 -12.281700 -11.21214 -1.128100 -6.280291
## 3 -6.651499 -14.82540 -3.552027 -6.325155
## 4 -12.473814 -13.54999 -2.420064 -6.289021
## 5 -11.711911 -13.73429 -2.595443 -2.278541
## 6 -10.442107 -14.43848 -3.406845 -2.307167
## log_scale_apap_Chan1997_1.4_0 log_scale_apap_glu_Chan1997_1.4_0
## 1 -0.34147420 0.3524540
## 2 0.87395974 -0.2202193
## 3 0.28419894 -1.2117128
## 4 -0.23704451 -0.5141464
## 5 -0.00256986 0.7640506
## 6 -1.39510256 -0.8312082
## log_scale_apap_sul_Chan1997_1.4_0 log_scale_apap_cys_Chan1997_1.4_0
## 1 -0.09054225 0.4909674
## 2 -2.05176756 0.1713548
## 3 1.59152206 0.9567482
## 4 -0.64296093 0.6822799
## 5 -0.28381384 1.2033634
## 6 -0.42447087 -0.9383802
## log_srel_apap_Chan1997_1.4_0 log_s0_apap_Chan1997_1.4_0
## 1 0.2569589 -1.9014794
## 2 -0.6194964 1.4224062
## 3 0.4884257 0.4620627
## 4 0.4431686 0.3437440
## 5 -1.0774869 -0.7111688
## 6 -0.1059288 -0.4287323
## log_srel_apap_glu_Chan1997_1.4_0 log_s0_apap_glu_Chan1997_1.4_0
## 1 -0.9569274 -0.1547545
## 2 0.6830946 -0.9358321
## 3 0.1252113 1.4079458
## 4 -1.0505252 -1.1969937
## 5 1.3007331 -1.5656455
## 6 2.1081213 -1.7762023
## log_srel_apap_sul_Chan1997_1.4_0 log_s0_apap_sul_Chan1997_1.4_0
## 1 -2.5583855 -1.1857697
## 2 0.2056214 -1.4048882
## 3 -0.5206815 -3.0952762
## 4 1.7018681 -2.4248608
## 5 -1.9290153 -0.8441127
## 6 0.6754411 0.6658790
## log_srel_apap_cys_Chan1997_1.4_0 log_s0_apap_cys_Chan1997_1.4_0
## 1 -0.80413329 -1.76384101
## 2 -1.03467092 -2.63912882
## 3 0.12948514 -0.06831229
## 4 -1.23494485 -0.35465510
## 5 -1.45640160 -2.21684171
## 6 -0.08965077 -2.18967728
## log_scale_apap_Chiew2010_5.6_0 log_scale_apap_glu_Chiew2010_5.6_0
## 1 -0.73225919 0.29875836
## 2 0.07746059 -0.75638815
## 3 -0.41557105 0.66860572
## 4 -0.34446349 -0.09935669
## 5 -0.06790696 0.17262465
## 6 -0.36695303 0.98236576
## log_scale_apap_sul_Chiew2010_5.6_0 log_srel_apap_Chiew2010_5.6_0
## 1 -0.2751353 0.2001658
## 2 -0.8853731 -1.1763805
## 3 -1.4579817 1.4459862
## 4 1.9013839 -1.1825723
## 5 0.8957957 0.4329237
## 6 -0.5285034 -1.6273736
## log_s0_apap_Chiew2010_5.6_0 log_srel_apap_glu_Chiew2010_5.6_0
## 1 -1.3083801 0.2537679
## 2 -0.9910735 0.4781154
## 3 -1.3495770 -0.2456441
## 4 -0.4922068 -0.4925589
## 5 1.2760837 1.4551869
## 6 -1.4419611 -2.2427514
## log_s0_apap_glu_Chiew2010_5.6_0 log_srel_apap_sul_Chiew2010_5.6_0
## 1 -0.9917074 0.30899346
## 2 -0.1322426 0.20086298
## 3 -2.4500674 1.13602514
## 4 -2.1513778 0.97349583
## 5 0.3285409 0.09686059
## 6 -1.1718316 -0.32388955
## log_s0_apap_sul_Chiew2010_5.6_0 log_scale_apap_Critchley2005_1.4_0
## 1 -0.2317920 0.75661973
## 2 -0.1080672 0.40330247
## 3 -1.2538891 -1.84876130
## 4 0.5970107 0.09984706
## 5 0.3377952 -1.08319140
## 6 -0.9579787 0.15268520
## log_scale_apap_glu_Critchley2005_1.4_0
## 1 -1.6075093
## 2 -0.4203937
## 3 -0.4089714
## 4 0.6322839
## 5 0.5396460
## 6 0.6909803
## log_scale_apap_sul_Critchley2005_1.4_0
## 1 1.1200586
## 2 1.3772116
## 3 -0.8219418
## 4 -0.1005233
## 5 -1.1104974
## 6 -0.4086092
## log_scale_apap_cys_Critchley2005_1.4_0 log_srel_apap_Critchley2005_1.4_0
## 1 0.7914627 0.5229484
## 2 -0.2206032 1.7093085
## 3 -0.7028494 -0.4162756
## 4 -0.3442296 -1.5224635
## 5 -0.7764242 -1.5065604
## 6 -1.2219040 -0.4096500
## log_s0_apap_Critchley2005_1.4_0 log_srel_apap_glu_Critchley2005_1.4_0
## 1 2.8885995 -0.6958732
## 2 0.2769623 1.6113013
## 3 1.9668600 0.8816571
## 4 0.4696681 0.6190646
## 5 -0.1120432 0.9772918
## 6 2.0007139 -1.2027400
## log_s0_apap_glu_Critchley2005_1.4_0
## 1 -1.462525446
## 2 -0.876828265
## 3 -2.072184801
## 4 -0.374247260
## 5 -0.093214581
## 6 0.005710159
## log_srel_apap_sul_Critchley2005_1.4_0
## 1 0.09474047
## 2 1.05831194
## 3 1.37006018
## 4 0.49068959
## 5 -0.15389290
## 6 -0.21423115
## log_s0_apap_sul_Critchley2005_1.4_0
## 1 -1.5436022
## 2 -1.9787239
## 3 1.1611746
## 4 -0.5722336
## 5 -0.7604901
## 6 -0.5263634
## log_srel_apap_cys_Critchley2005_1.4_0
## 1 0.6346692
## 2 -0.9788842
## 3 -1.7831523
## 4 0.6156580
## 5 -1.3146436
## 6 0.4575737
## log_s0_apap_cys_Critchley2005_1.4_0 log_scale_apap_Rawlins1977_0_1
## 1 0.4977260 0.021621889
## 2 -0.6198171 0.021845011
## 3 -1.8375393 0.008699087
## 4 -1.2241609 0.022133087
## 5 -1.1091981 0.026012429
## 6 -0.1110733 0.024257471
## log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1 2.454559 -8.977730
## 2 2.454535 -8.977814
## 3 2.483257 -8.996088
## 4 2.453265 -8.974551
## 5 2.437694 -8.941497
## 6 2.442176 -8.947284
## log_scale_apap_Rawlins1977_1_0 log_srel_apap_Rawlins1977_1_0
## 1 0.1678273 2.582696
## 2 0.1692129 2.581302
## 3 0.1388991 2.595842
## 4 0.1692214 2.580949
## 5 0.2307545 2.589881
## 6 0.2203842 2.591367
## log_s0_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 -14.66563 0.4137086
## 2 -16.05141 0.4157914
## 3 -15.53460 0.3500338
## 4 -15.71205 0.4157741
## 5 -16.02509 0.4775457
## 6 -15.33039 0.4570333
## log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1 2.341414 -15.33500
## 2 2.341451 -16.63546
## 3 2.347062 -16.04883
## 4 2.341416 -16.30266
## 5 2.346485 -16.51624
## 6 2.347016 -15.82922
## log_scale_apap_Rawlins1977_0.5_0 log_srel_apap_Rawlins1977_0.5_0
## 1 -1.5900253 1.9345084
## 2 -0.2137423 0.2236760
## 3 0.8627395 0.6101951
## 4 1.1149818 -0.5099961
## 5 -0.4089437 1.0647446
## 6 -0.1155511 0.7124745
## log_s0_apap_Rawlins1977_0.5_0
## 1 -2.5199219
## 2 -0.2345078
## 3 -1.2819850
## 4 1.5677611
## 5 -0.7884726
## 6 -0.3613277
## logAPAPCYS_HLM_CL: 1.2793391679861e-06
## logAPAPCYS_Km: 0.0438798488903149
## logAPAPGLU_HLM_CL: 2.32752164063581
## logAPAPGLU_Km: 21.520668349057
## logAPAPSUL_HLM_CL: 2.06404786189464e-05
## logKa_apap: 0.00184517865281865
## log_s0_apap_Chan1997_1.4_0: 0.149347509305303
## log_s0_apap_Chiew2010_5.6_0: 0.270257502535883
## log_s0_apap_Critchley2005_1.4_0: 17.9681281235639
## log_s0_apap_cys_Chan1997_1.4_0: 0.171385306079452
## log_s0_apap_cys_Critchley2005_1.4_0: 1.64497626822607
## log_s0_apap_glu_Chan1997_1.4_0: 0.856625476079955
## log_s0_apap_glu_Chiew2010_5.6_0: 0.370942798641939
## log_s0_apap_glu_Critchley2005_1.4_0: 0.231650514577204
## log_s0_apap_Rawlins1977_0_1: 0.000126188936947135
## log_s0_apap_Rawlins1977_0.5_0: 0.0804658932475238
## log_s0_apap_Rawlins1977_1_0: 4.27362459933834e-07
## log_s0_apap_Rawlins1977_2_0: 2.18822833245219e-07
## log_s0_apap_sul_Chan1997_1.4_0: 0.30551094502221
## log_s0_apap_sul_Chiew2010_5.6_0: 0.793111039544634
## log_s0_apap_sul_Critchley2005_1.4_0: 0.21361025170362
## log_scale_apap_Chan1997_1.4_0: 0.710721805251837
## log_scale_apap_Chiew2010_5.6_0: 0.480821497303708
## log_scale_apap_Critchley2005_1.4_0: 2.13106047495752
## log_scale_apap_cys_Chan1997_1.4_0: 1.63389603961365
## log_scale_apap_cys_Critchley2005_1.4_0: 2.20662164565056
## log_scale_apap_glu_Chan1997_1.4_0: 1.42255421536611
## log_scale_apap_glu_Chiew2010_5.6_0: 1.34818380327456
## log_scale_apap_glu_Critchley2005_1.4_0: 0.200386104365974
## log_scale_apap_Rawlins1977_0_1: 1.02185733567333
## log_scale_apap_Rawlins1977_0.5_0: 0.203920455093923
## log_scale_apap_Rawlins1977_1_0: 1.18273235599321
## log_scale_apap_Rawlins1977_2_0: 1.51241638200544
## log_scale_apap_sul_Chan1997_1.4_0: 0.913435744877439
## log_scale_apap_sul_Chiew2010_5.6_0: 0.759469382839319
## log_scale_apap_sul_Critchley2005_1.4_0: 3.06503365611115
## log_srel_apap_Chan1997_1.4_0: 1.29299192947592
## log_srel_apap_Chiew2010_5.6_0: 1.22160522476788
## log_srel_apap_Critchley2005_1.4_0: 1.68699428030717
## log_srel_apap_cys_Chan1997_1.4_0: 0.447475592070443
## log_srel_apap_cys_Critchley2005_1.4_0: 1.88639794692755
## log_srel_apap_glu_Chan1997_1.4_0: 0.384071160069748
## log_srel_apap_glu_Chiew2010_5.6_0: 1.28887257124136
## log_srel_apap_glu_Critchley2005_1.4_0: 0.498638857365226
## log_srel_apap_Rawlins1977_0_1: 11.6412981287543
## log_srel_apap_Rawlins1977_0.5_0: 6.92064127305638
## log_srel_apap_Rawlins1977_1_0: 13.2327664425859
## log_srel_apap_Rawlins1977_2_0: 10.3959247377831
## log_srel_apap_sul_Chan1997_1.4_0: 0.0774296530505137
## log_srel_apap_sul_Chiew2010_5.6_0: 1.36205346050779
## log_srel_apap_sul_Critchley2005_1.4_0: 1.099373499422959 Introduce error model and scaling factors in the dynamic model - not so bad, but not good either
- some scaling factors look as if they could help
- dynamics is described quite ok
- in apap_sul, model is still too slow
- in apap_glu, model is a bit too slow
- chiew 2010, apap_glu is better than in try 10.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters9 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters9 <- pars[!(names(pars)%in%c(free_parameters9,names(f)[1]))] %>% names
mydatalist <- data %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables9 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters9 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters9 <- c(free_parameters9, scale_parameters9)
error_model9 <- c(apap = "srel_apap*apap^2 +s0_apap",
apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters9 <- setdiff(getSymbols(error_model9), names(error_model9)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters9] <- paste0("exp(log", free_parameters9, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters9, x = scale_parameters9, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters9)
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters9, x = error_parameters9, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters9)
errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p9 <- NULL
for(i in 1:length(p_list)) { p9 <<- p9 + p_list[[i]]}
g9 <- Y(observables9, x)#, parameters = c(free_parameters9, scale_parameters9))
err9 <- Y(error_model9, g9)
obj9 <- normL2(mydatalist, (g9*x*p9), errmodel = err9)
pouter9 <- rep(0, length(getParameters(obj9))) %>% set_names(getParameters(obj9))
pouter9[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# obj9(pouter9)
# job9 <- runbg({myfit <- mstrust(objfun = obj9, center = pouter9, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht2", filename = "job9", input = global_env_without(c("mypred", "job", "myfit")))
# save(job9, file = "job9.rda")
# job9$check()
# # load("job9.rda")
# myfit9 <- job9$get()$ruprecht2
# save(myfit9, file = "myfit9.rda")
# job9$purge()
load("myfit9.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit9 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit9 %>% as.parframe())#+scale_y_log10()
mypred9 <- (g9*x*p9)(mytimes*4, myfit9 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot9 <- plotCombined(mypred9, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot9)
# myplot9
# obj9(myfit9 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.})## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 1 -2477.716 FALSE 100 -5.099331 7.999288
## 2 35 -2477.716 FALSE 100 -5.032580 7.301881
## 3 49 -2477.715 FALSE 100 -5.087414 7.628175
## 4 179 -2477.715 FALSE 100 -5.087632 7.116519
## 5 64 -2477.715 TRUE 99 -5.073408 7.132834
## 6 127 -2477.715 FALSE 100 -5.094321 6.967860
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -4.797299 -17.83559 -6.250228 -7.073747
## 2 -4.797371 -17.08140 -6.250285 -7.073740
## 3 -4.797351 -17.45008 -6.250078 -7.073678
## 4 -4.797346 -16.93769 -6.250345 -7.073796
## 5 -4.797347 -16.94634 -6.250387 -7.073790
## 6 -4.797427 -16.79246 -6.250201 -7.073706
## log_scale_apap_Chan1997_14_0 log_scale_apap_glu_Chan1997_14_0
## 1 0.3201883 10.240255
## 2 0.3201823 9.476110
## 3 0.3201498 9.857194
## 4 0.3202139 9.345829
## 5 0.3202098 9.347900
## 6 0.3201586 9.203812
## log_scale_apap_sul_Chan1997_14_0 log_scale_apap_cys_Chan1997_14_0
## 1 -0.5703412 7.541992
## 2 -0.5702552 6.787785
## 3 -0.5703223 7.156501
## 4 -0.5702481 6.644072
## 5 -0.5702448 6.652709
## 6 -0.5702149 6.498853
## log_srel_apap_Chan1997_14_0 log_s0_apap_Chan1997_14_0
## 1 -1.066589 -8.553130
## 2 -1.066910 -8.552872
## 3 -1.064453 -8.553918
## 4 -1.069107 -8.552475
## 5 -1.070049 -8.552360
## 6 -1.066075 -8.552996
## log_srel_apap_glu_Chan1997_14_0 log_s0_apap_glu_Chan1997_14_0
## 1 3.937521 -7.368639
## 2 3.937853 -7.368713
## 3 3.937651 -7.368663
## 4 3.937780 -7.368695
## 5 3.937766 -7.368738
## 6 3.937715 -7.368682
## log_srel_apap_sul_Chan1997_14_0 log_s0_apap_sul_Chan1997_14_0
## 1 4.178282 -9.337702
## 2 4.178646 -9.337195
## 3 4.178415 -9.337887
## 4 4.178152 -9.337426
## 5 4.178408 -9.337071
## 6 4.178751 -9.337377
## log_srel_apap_cys_Chan1997_14_0 log_s0_apap_cys_Chan1997_14_0
## 1 -4.014309 -9.955108
## 2 -3.320616 -9.955241
## 3 -3.635824 -9.955119
## 4 -3.161390 -9.955246
## 5 -3.162397 -9.955206
## 6 -2.983036 -9.955250
## log_scale_apap_Chiew2010_56_0 log_scale_apap_glu_Chiew2010_56_0
## 1 0.1542391 10.561511
## 2 0.1542455 9.797452
## 3 0.1541823 10.178496
## 4 0.1542872 9.667135
## 5 0.1542769 9.669024
## 6 0.1542120 9.525162
## log_scale_apap_sul_Chiew2010_56_0 log_srel_apap_Chiew2010_56_0
## 1 -0.7426937 0.3778347
## 2 -0.7426260 0.3779957
## 3 -0.7426694 0.3771786
## 4 -0.7426272 0.3780723
## 5 -0.7426270 0.3777571
## 6 -0.7425874 0.3779534
## log_s0_apap_Chiew2010_56_0 log_srel_apap_glu_Chiew2010_56_0
## 1 -7.069136 2.369965
## 2 -7.067144 2.370166
## 3 -7.067803 2.370305
## 4 -7.066808 2.371076
## 5 -7.066785 2.371267
## 6 -7.067703 2.370746
## log_s0_apap_glu_Chiew2010_56_0 log_srel_apap_sul_Chiew2010_56_0
## 1 -5.863964 1.413941
## 2 -5.863894 1.413643
## 3 -5.864205 1.414305
## 4 -5.864437 1.413383
## 5 -5.865029 1.413175
## 6 -5.864204 1.413946
## log_s0_apap_sul_Chiew2010_56_0 log_scale_apap_Critchley2005_14_0
## 1 -6.186949 0.2345061
## 2 -6.187438 0.2344554
## 3 -6.187013 0.2343445
## 4 -6.187419 0.2344923
## 5 -6.187215 0.2345560
## 6 -6.187368 0.2344314
## log_scale_apap_glu_Critchley2005_14_0
## 1 10.314313
## 2 9.550148
## 3 9.931238
## 4 9.419870
## 5 9.421957
## 6 9.277840
## log_scale_apap_sul_Critchley2005_14_0
## 1 -1.148680
## 2 -1.148614
## 3 -1.148658
## 4 -1.148615
## 5 -1.148612
## 6 -1.148576
## log_scale_apap_cys_Critchley2005_14_0 log_srel_apap_Critchley2005_14_0
## 1 7.232789 2.288037
## 2 6.478583 2.284297
## 3 6.847288 2.279003
## 4 6.334869 2.284317
## 5 6.343503 2.285165
## 6 6.189646 2.284663
## log_s0_apap_Critchley2005_14_0 log_srel_apap_glu_Critchley2005_14_0
## 1 -8.080480 -5.922499
## 2 -8.074298 -5.172423
## 3 -8.068500 -5.541084
## 4 -8.073987 -5.031696
## 5 -8.075026 -5.034967
## 6 -8.075041 -4.887071
## log_s0_apap_glu_Critchley2005_14_0 log_srel_apap_sul_Critchley2005_14_0
## 1 -6.304482 4.463234
## 2 -6.304953 4.464450
## 3 -6.304470 4.464503
## 4 -6.304871 4.464421
## 5 -6.304808 4.464391
## 6 -6.304898 4.464320
## log_s0_apap_sul_Critchley2005_14_0 log_srel_apap_cys_Critchley2005_14_0
## 1 -7.946701 -2.811800
## 2 -7.945526 -2.085599
## 3 -7.945487 -2.428260
## 4 -7.945550 -1.991045
## 5 -7.945590 -1.941678
## 6 -7.945679 -1.917429
## log_s0_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1 -10.12659 -0.1015749
## 2 -10.12657 -0.1015765
## 3 -10.12661 -0.1016133
## 4 -10.12658 -0.1015472
## 5 -10.12661 -0.1015634
## 6 -10.12657 -0.1015950
## log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1 2.881313 -10.45464
## 2 2.880629 -10.44504
## 3 2.880721 -10.44162
## 4 2.880528 -10.44397
## 5 2.879162 -10.41924
## 6 2.880092 -10.43567
## log_scale_apap_Rawlins1977_1_0 log_srel_apap_Rawlins1977_1_0
## 1 -0.03506085 3.061212
## 2 -0.03507190 3.062111
## 3 -0.03511765 3.062435
## 4 -0.03503499 3.062118
## 5 -0.03503313 3.062624
## 6 -0.03510720 3.062807
## log_s0_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 -8.476342 0.03209758
## 2 -8.477932 0.03208771
## 3 -8.478216 0.03205411
## 4 -8.478040 0.03212816
## 5 -8.478879 0.03214222
## 6 -8.479337 0.03205324
## log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1 2.456148 -8.369891
## 2 2.457079 -8.371340
## 3 2.457015 -8.370964
## 4 2.457062 -8.371123
## 5 2.457133 -8.372795
## 6 2.456935 -8.370860
## log_scale_apap_Rawlins1977_05_0 log_srel_apap_Rawlins1977_05_0
## 1 -0.2098564 -7.088126
## 2 -0.2098571 -6.338162
## 3 -0.2098887 -6.706408
## 4 -0.2098285 -6.195749
## 5 -0.2098314 -6.201175
## 6 -0.2098756 -6.053133
## log_s0_apap_Rawlins1977_05_0
## 1 -8.283322
## 2 -8.284305
## 3 -8.283902
## 4 -8.284432
## 5 -8.284290
## 6 -8.284218
## logAPAPCYS_HLM_CL: 1.795146994404e-08
## logAPAPCYS_Km: 0.00193001355192211
## logAPAPGLU_HLM_CL: 0.00610082435223236
## logAPAPGLU_Km: 2978.83487681056
## logAPAPSUL_HLM_CL: 0.00825200823079057
## logKa_apap: 0.000847053225794413
## log_s0_apap_Chan1997_14_0: 0.00019294018495771
## log_s0_apap_Chiew2010_56_0: 0.00085096838208337
## log_s0_apap_Critchley2005_14_0: 0.000309522401143699
## log_s0_apap_cys_Chan1997_14_0: 4.74844786235421e-05
## log_s0_apap_cys_Critchley2005_14_0: 4.00018130602822e-05
## log_s0_apap_glu_Chan1997_14_0: 0.000630726272006818
## log_s0_apap_glu_Chiew2010_56_0: 0.00283996453762311
## log_s0_apap_glu_Critchley2005_14_0: 0.00182809286937598
## log_s0_apap_Rawlins1977_0_1: 2.8814215684246e-05
## log_s0_apap_Rawlins1977_05_0: 0.000252696390036355
## log_s0_apap_Rawlins1977_1_0: 0.000208339338154456
## log_s0_apap_Rawlins1977_2_0: 0.000231740837081231
## log_s0_apap_sul_Chan1997_14_0: 8.80414871966006e-05
## log_s0_apap_sul_Chiew2010_56_0: 0.00205609044641647
## log_s0_apap_sul_Critchley2005_14_0: 0.000353827456947821
## log_scale_apap_Chan1997_14_0: 1.37738706359055
## log_scale_apap_Chiew2010_56_0: 1.16676983632714
## log_scale_apap_Critchley2005_14_0: 1.26428417568878
## log_scale_apap_cys_Chan1997_14_0: 1885.5826917627
## log_scale_apap_cys_Critchley2005_14_0: 1384.0773914102
## log_scale_apap_glu_Chan1997_14_0: 28008.2659198081
## log_scale_apap_glu_Chiew2010_56_0: 38619.4394146001
## log_scale_apap_glu_Critchley2005_14_0: 30161.2433788063
## log_scale_apap_Rawlins1977_0_1: 0.903413505390761
## log_scale_apap_Rawlins1977_05_0: 0.81070064561698
## log_scale_apap_Rawlins1977_1_0: 0.965546657143368
## log_scale_apap_Rawlins1977_2_0: 1.03261825838704
## log_scale_apap_sul_Chan1997_14_0: 0.565332533789827
## log_scale_apap_sul_Chiew2010_56_0: 0.475830435424811
## log_scale_apap_sul_Critchley2005_14_0: 0.317054858672663
## log_srel_apap_Chan1997_14_0: 0.344180473770218
## log_srel_apap_Chiew2010_56_0: 1.45912166820889
## log_srel_apap_Critchley2005_14_0: 9.85557305660906
## log_srel_apap_cys_Chan1997_14_0: 0.0180554247494083
## log_srel_apap_cys_Critchley2005_14_0: 0.0600966934703123
## log_srel_apap_glu_Chan1997_14_0: 51.2912931730463
## log_srel_apap_glu_Chiew2010_56_0: 10.6970164613051
## log_srel_apap_glu_Critchley2005_14_0: 0.00267849706185607
## log_srel_apap_Rawlins1977_0_1: 17.8376851856982
## log_srel_apap_Rawlins1977_05_0: 0.000834960350013315
## log_srel_apap_Rawlins1977_1_0: 21.3534269462142
## log_srel_apap_Rawlins1977_2_0: 11.6598121177031
## log_srel_apap_sul_Chan1997_14_0: 65.2536831006089
## log_srel_apap_sul_Chiew2010_56_0: 4.11213078704973
## log_srel_apap_sul_Critchley2005_14_0: 86.767695781153310 Introduce error model only - quite nice, problems in chiew2100 and in apap_sul
- some scaling factors look as if they could help
- dynamics is described quite ok
- in apap_sul, model is still too slow
- in apap_glu, model is a bit too slow
- chiew 2010, apap_sul is better than in tr9
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters10 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters10 <- pars[!(names(pars)%in%c(free_parameters10,names(f)[1]))] %>% names
mydatalist <- data %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables10 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters10 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters10 <- c(free_parameters10, scale_parameters10)
error_model10 <- c(apap = "srel_apap*apap^2 +s0_apap",
apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
error_parameters10 <- setdiff(getSymbols(error_model10), names(error_model10)) %>% set_names(.,.)
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters10] <- paste0("exp(log", free_parameters10, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters10, x = scale_parameters10, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters10)
scales[1:length(scales)] <- "1"
errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters10, x = error_parameters10, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters10)
errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p10 <- NULL
for(i in 1:length(p_list)) { p10 <<- p10 + p_list[[i]]}
g10 <- Y(observables10, x)#, parameters = c(free_parameters10, scale_parameters10))
err10 <- Y(error_model10, g10)
obj10 <- normL2(mydatalist, (g10*x*p10), errmodel = err10)
pouter10 <- rep(0, length(getParameters(obj10))) %>% set_names(getParameters(obj10))
pouter10[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# job10 <- runbg({myfit <- mstrust(objfun = obj10, center = pouter10, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht1", filename = "job10", input = global_env_without(c("mypred", "job", "myfit")))
# save(job10, file = "job10.rda")
# job10$check()# load("job10.rda")
# myfit10 <- job10$get()$ruprecht1
# save(myfit10, file = "myfit10.rda")
# job10$purge()
load("myfit10.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit10 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit10 %>% as.parframe())#+scale_y_log10()
mypred10 <- (g10*x*p10)(mytimes*4, myfit10 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot10 <- plotCombined(mypred10, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot10)
# myplot10## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 72 -2356.519 TRUE 99 -5.984038 -3.407580
## 2 23 -2356.519 TRUE 100 -5.984055 -3.407593
## 3 90 -2356.519 TRUE 98 -5.984039 -3.407588
## 4 144 -2356.519 FALSE 100 -5.984040 -3.407582
## 5 63 -2356.519 TRUE 78 -5.983866 -3.407695
## 6 135 -2356.519 FALSE 100 -5.984081 -3.407616
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -5.782424 -10.57408 -6.696947 -7.013798
## 2 -5.782408 -10.57408 -6.696946 -7.013783
## 3 -5.782429 -10.57408 -6.696951 -7.013733
## 4 -5.782424 -10.57408 -6.696946 -7.013786
## 5 -5.782736 -10.57412 -6.697111 -7.013753
## 6 -5.782389 -10.57409 -6.696947 -7.013793
## log_srel_apap_Chan1997_14_0 log_s0_apap_Chan1997_14_0
## 1 3.289950 -7.052483
## 2 3.289934 -7.052496
## 3 3.289921 -7.052506
## 4 3.289937 -7.052506
## 5 3.289810 -7.052222
## 6 3.289991 -7.052427
## log_srel_apap_glu_Chan1997_14_0 log_s0_apap_glu_Chan1997_14_0
## 1 4.404972 -7.517436
## 2 4.405154 -7.517363
## 3 4.405028 -7.517374
## 4 4.405060 -7.517341
## 5 4.403288 -7.516821
## 6 4.405348 -7.517387
## log_srel_apap_sul_Chan1997_14_0 log_s0_apap_sul_Chan1997_14_0
## 1 4.959255 -8.122914
## 2 4.958957 -8.123157
## 3 4.959125 -8.123015
## 4 4.959096 -8.123068
## 5 4.960016 -8.122510
## 6 4.958673 -8.123290
## log_srel_apap_cys_Chan1997_14_0 log_s0_apap_cys_Chan1997_14_0
## 1 -6.968630 -9.728692
## 2 -3.817880 -9.728688
## 3 -5.081099 -9.728682
## 4 -3.199117 -9.728689
## 5 -2.727214 -9.728656
## 6 -2.482580 -9.728732
## log_srel_apap_Chiew2010_56_0 log_s0_apap_Chiew2010_56_0
## 1 0.1017979 -7.141670
## 2 0.1020375 -7.141873
## 3 0.1016049 -7.141541
## 4 0.1018871 -7.141784
## 5 0.1016913 -7.141588
## 6 0.1019149 -7.141807
## log_srel_apap_glu_Chiew2010_56_0 log_s0_apap_glu_Chiew2010_56_0
## 1 5.497352 -6.479678
## 2 5.497414 -6.479689
## 3 5.497194 -6.479812
## 4 5.497293 -6.479720
## 5 5.496531 -6.479613
## 6 5.498219 -6.479841
## log_srel_apap_sul_Chiew2010_56_0 log_s0_apap_sul_Chiew2010_56_0
## 1 2.464006 -6.105482
## 2 2.463693 -6.105647
## 3 2.463902 -6.105591
## 4 2.463910 -6.105581
## 5 2.468499 -6.106767
## 6 2.463066 -6.105751
## log_srel_apap_Critchley2005_14_0 log_s0_apap_Critchley2005_14_0
## 1 2.350997 -6.825250
## 2 2.350973 -6.825248
## 3 2.350992 -6.825255
## 4 2.350957 -6.825260
## 5 2.350889 -6.824922
## 6 2.350855 -6.825138
## log_srel_apap_glu_Critchley2005_14_0 log_s0_apap_glu_Critchley2005_14_0
## 1 3.602125 -6.967212
## 2 3.602291 -6.967185
## 3 3.602170 -6.967135
## 4 3.602175 -6.967151
## 5 3.601192 -6.966273
## 6 3.602554 -6.967230
## log_srel_apap_sul_Critchley2005_14_0 log_s0_apap_sul_Critchley2005_14_0
## 1 4.017564 -7.716869
## 2 4.017535 -7.716892
## 3 4.017606 -7.716842
## 4 4.017561 -7.716888
## 5 4.017116 -7.717051
## 6 4.017688 -7.716863
## log_srel_apap_cys_Critchley2005_14_0 log_s0_apap_cys_Critchley2005_14_0
## 1 -8.051930 -10.10211
## 2 -4.902004 -10.10212
## 3 -6.164908 -10.10211
## 4 -4.283094 -10.10212
## 5 -3.755340 -10.10221
## 6 -3.566113 -10.10215
## log_srel_apap_Rawlins1977_0_1 log_s0_apap_Rawlins1977_0_1
## 1 2.541371 -8.725966
## 2 2.541383 -8.725981
## 3 2.541391 -8.725948
## 4 2.541383 -8.725955
## 5 2.541252 -8.725793
## 6 2.541297 -8.725894
## log_srel_apap_Rawlins1977_1_0 log_s0_apap_Rawlins1977_1_0
## 1 2.772759 -8.547667
## 2 2.772043 -8.546641
## 3 2.770963 -8.545399
## 4 2.772593 -8.547197
## 5 2.772166 -8.547219
## 6 2.772761 -8.547107
## log_srel_apap_Rawlins1977_2_0 log_s0_apap_Rawlins1977_2_0
## 1 2.485538 -8.345503
## 2 2.485019 -8.343667
## 3 2.485324 -8.344722
## 4 2.485179 -8.344683
## 5 2.485390 -8.344672
## 6 2.485404 -8.344505
## log_srel_apap_Rawlins1977_05_0 log_s0_apap_Rawlins1977_05_0
## 1 3.075085 -8.175010
## 2 3.075391 -8.175087
## 3 3.075176 -8.175017
## 4 3.074992 -8.174969
## 5 3.074871 -8.175119
## 6 3.074939 -8.175055
## logAPAPCYS_HLM_CL: 2.5570243395892e-05
## logAPAPCYS_Km: 0.00123467580981531
## logAPAPGLU_HLM_CL: 0.00251863572465911
## logAPAPGLU_Km: 0.033121273391679
## logAPAPSUL_HLM_CL: 0.00308123765707667
## logKa_apap: 0.000899386426280096
## log_s0_apap_Chan1997_14_0: 0.000865258189660525
## log_s0_apap_Chiew2010_56_0: 0.000791429248666158
## log_s0_apap_Critchley2005_14_0: 0.00108600472415769
## log_s0_apap_cys_Chan1997_14_0: 5.9550106934986e-05
## log_s0_apap_cys_Critchley2005_14_0: 4.09928708157083e-05
## log_s0_apap_glu_Chan1997_14_0: 0.000543524495738245
## log_s0_apap_glu_Chiew2010_56_0: 0.00153430475174298
## log_s0_apap_glu_Critchley2005_14_0: 0.000942275918358174
## log_s0_apap_Rawlins1977_0_1: 0.000162315884677548
## log_s0_apap_Rawlins1977_05_0: 0.000281603560301461
## log_s0_apap_Rawlins1977_1_0: 0.000193997116676406
## log_s0_apap_Rawlins1977_2_0: 0.000237462036719598
## log_s0_apap_sul_Chan1997_14_0: 0.000296662808440626
## log_s0_apap_sul_Chiew2010_56_0: 0.00223060565323054
## log_s0_apap_sul_Critchley2005_14_0: 0.000445252327814248
## log_srel_apap_Chan1997_14_0: 26.8415245012767
## log_srel_apap_Chiew2010_56_0: 1.10715971018415
## log_srel_apap_Critchley2005_14_0: 10.4960243104907
## log_srel_apap_cys_Chan1997_14_0: 0.00094094138090811
## log_srel_apap_cys_Critchley2005_14_0: 0.000318486578636513
## log_srel_apap_glu_Chan1997_14_0: 81.8568173051488
## log_srel_apap_glu_Chiew2010_56_0: 244.044756600744
## log_srel_apap_glu_Critchley2005_14_0: 36.6760946617877
## log_srel_apap_Rawlins1977_0_1: 12.6970698875144
## log_srel_apap_Rawlins1977_05_0: 21.6517119661049
## log_srel_apap_Rawlins1977_1_0: 16.0027176603338
## log_srel_apap_Rawlins1977_2_0: 12.0075810603921
## log_srel_apap_sul_Chan1997_14_0: 142.487667782381
## log_srel_apap_sul_Chiew2010_56_0: 11.7517960003797
## log_srel_apap_sul_Critchley2005_14_0: 55.5656002676397Data revisited
Some sigmas of the data are NA, try to recover an estimate for the sigmas from the other sigmas.
Take out some outliers for the fitting. Those are: study, name, time, reason 1. Chan1997, apap_cys, 5400, don’t know what went wrong with this one, but it just doesn’t fit reasonably in the time course 3. Chan1997, apap, 18000, The sigma is a few orders of magnitude lower. Mirjam said it might be that in this point less people were measured (eg 2) and they had nearly the same value. Then of course, sigma would be very smal
# data <- data %>%
# filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
# filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
# {.}
data %>% select(-n) %>% filter(study %>% str_detect("Chiew"))%>% as.datalist() %>% plotData()
data_with_errors <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name"), blather = T)
data_with_errors %>% select(study, name, D_apap, s0, srel) %>% unique()## study name D_apap s0 srel
## 1 Chan1997 apap 1.4 -14.395145 -4.516431
## 13 Chan1997 apap_sul 1.4 -15.377176 -2.609451
## 26 Chan1997 apap_glu 1.4 -14.970003 -1.709818
## 39 Chan1997 apap_cys 1.4 -19.186945 -13.879472
## 51 Chiew2010 apap 5.6 -24.503629 -1.232431
## 65 Chiew2010 apap_glu 5.6 -9.710964 -2.303673
## 79 Chiew2010 apap_sul 5.6 -10.576758 -4.220755
## 93 Critchley2005 apap 1.4 -11.298969 -16.644666
## 108 Critchley2005 apap_sul 1.4 -16.328137 -2.802784
## 123 Critchley2005 apap_glu 1.4 -15.826025 -2.457868
## 138 Critchley2005 apap_cys 1.4 -18.851664 -18.687321
## 153 Rawlins1977 apap 0 -13.826728 -2.924789
## 163 Rawlins1977 apap 1 -28.577775 -1.968079
## 171 Rawlins1977 apap 2 -27.097492 -1.002046
## 179 Rawlins1977 apap 0.5 -27.374147 -2.115781What’s the problem with the Chiew-dataset? The error model doesn’t work out, somehow
data_with_errors %>%
filter(study %>% str_detect("Chiew"), name %in% "apap") %>%
ggplot(aes(x=value)) +
geom_point(aes(y=sigmaLS^2*(n), color = log(time))) +
geom_line(aes(y=sigma^2*n)) +
geom_ribbon(aes(ymin=cbLower95, ymax=cbUpper95), alpha=.3) +
geom_ribbon(aes(ymin=cbLower68, ymax=cbUpper68), alpha=.3) +
ylab("variance") +
facet_wrap(~condidnt, scales = "free") +
scale_y_log10() +
theme_dMod() +
scale_color_continuous( low = "#98f5ff", high = "#4c4cdb")myplot <- data_with_errors %>%
rename(sigma_fitted = sigma) %>%
left_join(data) %>%
gather("which_sig", "sigma", sigma, sigma_fitted) %>%
mutate(condition = paste0(study, D_apap)) %>%
ggplot(aes(x = time,y = log10(sigma))) +
geom_line(aes(color = condition, linetype = which_sig)) +
# geom_point(aes(color = condition, shape = which_sig))+
facet_wrap("name", scales = "free") ## Joining, by = c("study", "n", "time", "name", "value", "D_apap", "Ave_apap")myplot %>% plotly::ggplotly()
# myplotFit error model first, then fit the model with free scaling factors
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters_e1 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters_e1 <- pars[!(names(pars)%in%c(free_parameters_e1,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_e1 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_e1 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_e1 <- c(free_parameters_e1, scale_parameters_e1)
# error_model_e1 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_e1 <- setdiff(getSymbols(error_model_e1), names(error_model_e1)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_e1] <- paste0("exp(log", free_parameters_e1, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_e1, x = scale_parameters_e1, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_e1)
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_e1, x = error_parameters_e1, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_e1)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_e1 <- NULL
for(i in 1:length(p_list)) { p_e1 <<- p_e1 + p_list[[i]]}
g_e1 <- Y(observables_e1, x)#, parameters = c(free_parameters_e1, scale_parameters_e1))
# err_e1 <- Y(error_model_e1, g_e1)
obj_e1 <- normL2(mydatalist, (g_e1*x*p_e1))#, errmodel = err_e1)
pouter_e1 <- rep(0, length(getParameters(obj_e1))) %>% set_names(getParameters(obj_e1))
# pouter_e1[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
#job_e1 <- runbg({myfit <- mstrust(objfun = obj_e1, center = pouter_e1, sd = 3, studyname = "methacetin", cores = 24, fits = 200); myfit}, machine = "ruprecht1", filename = "job_e1", input = global_env_without(c("mypred", "job", "myfit")))
# save(job_e1, file = "job_e1.rda")
# job_e1$check()# load("job_e1_1_result.RData")
# myfit_e1 <- .runbgOutput
# load("job_e1.rda")
# myfit_e1 <- job_e1$get()
# save(myfit_e1, file = "myfit_e1.rda")
# job_e1$purge()
load("myfit_e1.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit_e1 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_e1 %>% as.parframe())#+scale_y_log10()
mypred_e1 <- (g_e1*x*p_e1)(mytimes*4, myfit_e1 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot_e1 <- plotCombined(mypred_e1, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_e1)
# myplot_e1## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 135 1151.996 TRUE 69 -18.255493 -2.247950
## 2 70 1151.996 FALSE 100 -25.864555 -2.246878
## 3 102 1151.996 TRUE 87 -20.902541 -2.247412
## 4 171 1151.996 TRUE 61 -15.814734 -2.245022
## 5 191 1152.002 TRUE 60 -7.374898 -2.255883
## 6 81 1152.021 FALSE 100 -5.978911 -2.280757
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -4.439263 -9.981200 -7.521667 -7.664488
## 2 -4.439394 -9.980696 -7.521636 -7.664391
## 3 -4.439279 -9.981359 -7.521708 -7.664394
## 4 -4.439563 -9.979954 -7.521437 -7.664217
## 5 -4.482247 -9.973992 -7.522330 -7.663926
## 6 -4.629608 -9.955660 -7.524306 -7.662562
## log_scale_apap_Chan1997_14_0 log_scale_apap_glu_Chan1997_14_0
## 1 0.6599173 12.8364130
## 2 0.6598512 20.4463668
## 3 0.6599152 15.4839438
## 4 0.6597140 10.3980699
## 5 0.6598690 1.9483154
## 6 0.6599017 0.5291009
## log_scale_apap_sul_Chan1997_14_0 log_scale_apap_cys_Chan1997_14_0
## 1 -1.0771287 -0.4991120
## 2 -1.0770584 -0.4996198
## 3 -1.0770897 -0.4989502
## 4 -1.0770215 -0.5003547
## 5 -1.0342097 -0.5064348
## 6 -0.8868909 -0.5250863
## log_scale_apap_Chiew2010_56_0 log_scale_apap_glu_Chiew2010_56_0
## 1 0.4747969 14.360289
## 2 0.4746995 21.970031
## 3 0.4747278 17.007724
## 4 0.4745103 11.921363
## 5 0.4764708 3.475552
## 6 0.4813080 2.066233
## log_scale_apap_sul_Chiew2010_56_0 log_scale_apap_Critchley2005_14_0
## 1 -0.8302279 0.4550537
## 2 -0.8301517 0.4549816
## 3 -0.8301895 0.4550346
## 4 -0.8301000 0.4548393
## 5 -0.7857834 0.4549590
## 6 -0.6342932 0.4548708
## log_scale_apap_glu_Critchley2005_14_0
## 1 13.2497621
## 2 20.8597371
## 3 15.8973054
## 4 10.8114742
## 5 2.3620052
## 6 0.9437402
## log_scale_apap_sul_Critchley2005_14_0
## 1 -1.373527
## 2 -1.373459
## 3 -1.373493
## 4 -1.373425
## 5 -1.330606
## 6 -1.183280
## log_scale_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1 -0.7385315 0.1580270
## 2 -0.7390379 0.1579934
## 3 -0.7383715 0.1580312
## 4 -0.7397729 0.1579200
## 5 -0.7458288 0.1579958
## 6 -0.7644113 0.1579937
## log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 0.3682622 0.1562957
## 2 0.3682065 0.1562248
## 3 0.3682818 0.1562731
## 4 0.3680777 0.1560860
## 5 0.3680259 0.1564571
## 6 0.3675088 0.1570966
## log_scale_apap_Rawlins1977_05_0
## 1 0.1238156
## 2 0.1237872
## 3 0.1238273
## 4 0.1236985
## 5 0.1235803
## 6 0.1229841
## logAPAPCYS_HLM_CL: 4.6261539108133e-05
## logAPAPCYS_Km: 0.000541229613400164
## logAPAPGLU_HLM_CL: 1.17961442872092e-08
## logAPAPGLU_Km: 0.105615516447963
## logAPAPSUL_HLM_CL: 0.0118046378599577
## logKa_apap: 0.000469196966461811
## log_scale_apap_Chan1997_14_0: 1.93463224874753
## log_scale_apap_Chiew2010_56_0: 1.60768769704826
## log_scale_apap_Critchley2005_14_0: 1.57625801184958
## log_scale_apap_cys_Chan1997_14_0: 0.607069478303199
## log_scale_apap_cys_Critchley2005_14_0: 0.47781508039879
## log_scale_apap_glu_Chan1997_14_0: 375649.964718959
## log_scale_apap_glu_Chiew2010_56_0: 1724226.50632285
## log_scale_apap_glu_Critchley2005_14_0: 567934.927599969
## log_scale_apap_Rawlins1977_0_1: 1.1711977984144
## log_scale_apap_Rawlins1977_05_0: 1.13180719602812
## log_scale_apap_Rawlins1977_1_0: 1.44522093867267
## log_scale_apap_Rawlins1977_2_0: 1.16917192088704
## log_scale_apap_sul_Chan1997_14_0: 0.340572019805132
## log_scale_apap_sul_Chiew2010_56_0: 0.435949924368047
## log_scale_apap_sul_Critchley2005_14_0: 0.253212324969791Fit error model to data, then fit model without scaling
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters_e2 <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap"#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters_e2 <- pars[!(names(pars)%in%c(free_parameters_e2,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_e2 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_e2 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_e2 <- c(free_parameters_e2, scale_parameters_e2)
# error_model_e2 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_e2 <- setdiff(getSymbols(error_model_e2), names(error_model_e2)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_e2] <- paste0("exp(log", free_parameters_e2, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_e2, x = scale_parameters_e2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_e2)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_e2, x = error_parameters_e2, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_e2)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_e2 <- NULL
for(i in 1:length(p_list)) { p_e2 <<- p_e2 + p_list[[i]]}
g_e2 <- Y(observables_e2, x)#, parameters = c(free_parameters_e2, scale_parameters_e2))
# err_e2 <- Y(error_model_e2, g_e2)
obj_e2 <- normL2(mydatalist, (g_e2*x*p_e2))#, errmodel = err_e2)
pouter_e2 <- rep(0, length(getParameters(obj_e2))) %>% set_names(getParameters(obj_e2))
# pouter_e2[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# job_e2 <- runbg({myfit <- mstrust(objfun = obj_e2, center = pouter_e2, sd = 5, studyname = "methacetin", cores = 20, fits = 400); myfit}, machine = "ruprecht2", filename = "job_e2", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_e2, file = "job_e2.rda")
# job_e2$check()# load("job_e2_1_result.RData")
# myfit_e2 <- .runbgOutput
# load("job_e1.rda")
# myfit_e1 <- job_e1$get()
# save(myfit_e2, file = "myfit_e2.rda")
# job_e1$purge()
load("myfit_e2.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit_e2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_e2 %>% as.parframe())#+scale_y_log10()
mypred_e2 <- (g_e2*x*p_e2)(mytimes*4, myfit_e2 %>% as.parframe() %>% as.parvec %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
myplot_e2 <- plotCombined(mypred_e2, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_e2)
# myplot_e2## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 379 2402.825 TRUE 56 -6.195743 -3.728993
## 2 125 2402.825 TRUE 46 -6.195869 -3.729162
## 3 147 2402.826 TRUE 68 -6.195203 -3.728064
## 4 224 2402.826 TRUE 54 -6.195879 -3.729148
## 5 374 2402.826 TRUE 39 -6.195582 -3.728692
## 6 28 2402.826 TRUE 38 -6.195853 -3.729130
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logKa_apap
## 1 -5.922382 -10.65531 -7.895205 -7.156523
## 2 -5.922387 -10.65532 -7.895313 -7.156831
## 3 -5.922448 -10.65529 -7.894827 -7.156478
## 4 -5.922413 -10.65534 -7.895534 -7.156525
## 5 -5.922399 -10.65530 -7.895033 -7.156939
## 6 -5.922391 -10.65534 -7.895523 -7.156850
## logAPAPCYS_HLM_CL: 2.3575301587201e-05
## logAPAPCYS_Km: 0.000372525594735151
## logAPAPGLU_HLM_CL: 0.00203808917141075
## logAPAPGLU_Km: 0.0240169974041635
## logAPAPSUL_HLM_CL: 0.00267881249032963
## logKa_apap: 0.000779761067836627Tackle too slow dynamics of apap_*
The dynamics of e.g. apap_glu are too slow. Test two hypotheses:
- Increasing CLrenal_apap* fixes it
- multiplying a rate constant to Ave_apap_ -> A_lu_apap_ fixes it.
Compile a model with less sensitivities for faster integration
This model also comprises some rate constants to make some reactions go faster
# Multiply rates from "ve" to "lu" with additional free parameters
reactions_to_lung <- reactions
reactions_to_lung[reactions_to_lung$from %>% str_detect("(1*Ave_apap)"),3] %<>% paste0(" * (k_to_lung_", c("sul", "", "cys", "glu"),")")
reactions_to_lung
el_to_lung <- eqnlist()
for(i in 1:nrow(reactions_to_lung)) el_to_lung <- addReaction(el_to_lung, reactions_to_lung$from[i], reactions_to_lung$to[i], reactions_to_lung$rate[i], reactions_to_lung$description[i])
# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f_to_lung <- el_to_lung %>% as.eqnvec()
parameters <- all_pars[getParameters(x)] %>% names()
free_parameters_to_lung <- parameters[parameters %>% sapply(. %>% str_detect(c("CL", "Km", "Ka_apap$")) %>% any) ]
# myodemodel_to_lung <- odemodel(f_to_lung, modelname = "methacetin_to_lung", fixed = setdiff(parameters, c(free_parameters_to_lung, "Ave_apap"))) # Include "Ave_apap" because of a bug: You need the sensitivities to at least one initial value
save(myodemodel_to_lung, file = "methacetin_to_lung.rda")Fit the model with CLrenal*-parameters The hypothesis was that the slowly decaying dynamics of _glu, _cys and _sul could be fixed by removing them faster from the system due to clearance by the kidney. This didn’t work.
load("methacetin_to_lung.rda")
x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)
# get the only the parameters needed for x
pars_to_lung <- rep(1,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
pars <- c(all_pars, pars_to_lung)[getParameters(x_to_lung)]
free_parameters_CLrenal <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap",
"CLrenal_apap", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_apap_sul"
)
fixed_parameters_CLrenal <- pars[!(names(pars)%in%c(free_parameters_CLrenal,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_CLrenal <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_CLrenal <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_CLrenal <- c(free_parameters_CLrenal, scale_parameters_CLrenal)
# error_model_CLrenal <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_CLrenal <- setdiff(getSymbols(error_model_CLrenal), names(error_model_CLrenal)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_CLrenal] <- paste0("exp(log", free_parameters_CLrenal, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_CLrenal, x = scale_parameters_CLrenal, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_CLrenal)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_CLrenal, x = error_parameters_CLrenal, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_CLrenal)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_CLrenal <- NULL
for(i in 1:length(p_list)) { p_CLrenal <<- p_CLrenal + p_list[[i]]}
g_CLrenal <- Y(observables_CLrenal, x_to_lung)#, parameters = c(free_parameters_CLrenal, scale_parameters_CLrenal))
# err_CLrenal <- Y(error_model_CLrenal, g_CLrenal)
obj_CLrenal <- normL2(mydatalist, (g_CLrenal*x_to_lung*p_CLrenal))#, errmodel = err_CLrenal)
pouter_CLrenal <- rep(0, length(getParameters(obj_CLrenal))) %>% set_names(getParameters(obj_CLrenal))
pouter_CLrenal[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
pouter_CLrenal[names(pouter_CLrenal) %>% str_detect("CLr")] <- -4
pouter_CLrenal[1:2] <- c(1.5,5)
# obj_CLrenal(pouter_CLrenal)
# job_CLrenal <- runbg({myfit <- mstrust(objfun = obj_CLrenal, center = pouter_CLrenal, sd = 3, studyname = "methacetin", cores = 16, fits = 100, iterlim = 150); myfit}, machine = "ruprecht2", filename = "job_CLrenal", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_CLrenal, file = "job_CLrenal.rda")
# job_CLrenal$check()# myfit_CLrenal <- job_CLrenal$get()$ruprecht2
# # save(myfit_CLrenal, file = "myfit_CLrenal.rda")
# # job_CLrenal$purge()
load("myfit_CLrenal.rda")
loadDLL(x_to_lung)## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.somyfit_CLrenal %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_CLrenal %>% as.parframe())+scale_y_log10()
mypouter <- myfit_CLrenal %>% as.parframe() %>% as.parvec
# mypouter[1:2] <- mypouter[1:2]+1
mypred_CLrenal <- (g_CLrenal*x_to_lung*p_CLrenal)(mytimes*4, mypouter, deriv = F)
# myplot_CLrenal <- plotCombined(mypred_CLrenal, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_CLrenal <- plotCombined(mypred_CLrenal, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_CLrenal)
# myplot_CLrenal## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 22 1679.527 TRUE 55 -6.532773 -4.069886
## 2 7 1679.527 TRUE 46 -6.533129 -4.070574
## 3 60 1679.527 TRUE 59 -6.532886 -4.069996
## 4 31 1679.527 TRUE 60 -6.532866 -4.069961
## 5 70 1679.527 TRUE 55 -6.532908 -4.070019
## 6 71 1679.527 TRUE 66 -6.532916 -4.070109
## logAPAPSUL_HLM_CL logAPAPCYS_HLM_CL logAPAPCYS_Km logCLrenal_apap_sul
## 1 -5.629237 -10.64046 -6.378762 -5.099622
## 2 -5.629241 -10.64043 -6.374378 -5.099647
## 3 -5.629227 -10.64037 -6.382149 -5.099616
## 4 -5.629244 -10.64037 -6.383058 -5.099634
## 5 -5.629232 -10.64037 -6.383474 -5.099626
## 6 -5.629243 -10.64041 -6.380751 -5.099644
## logKa_apap logCLrenal_apap logCLrenal_apap_cys logCLrenal_apap_glu
## 1 -7.220531 -32.01661 -5.271552 -5.448181
## 2 -7.220536 -28.60979 -5.273667 -5.448045
## 3 -7.220632 -26.47175 -5.269753 -5.448253
## 4 -7.220556 -33.60964 -5.269310 -5.448251
## 5 -7.220625 -19.92413 -5.269103 -5.448267
## 6 -7.220526 -27.18079 -5.270518 -5.448182
## logAPAPCYS_HLM_CL: 2.39279181355532e-05
## logAPAPCYS_Km: 0.0016972225659718
## logAPAPGLU_HLM_CL: 0.00145496514549494
## logAPAPGLU_Km: 0.0170793409249032
## logAPAPSUL_HLM_CL: 0.00359131615388076
## logCLrenal_apap: 1.24556095423351e-14
## logCLrenal_apap_cys: 0.00513563209480641
## logCLrenal_apap_glu: 0.0043041255170805
## logCLrenal_apap_sul: 0.00609905356785048
## logKa_apap: 0.000731413761365896Allow rate from venous to lung to be faster - 2 local optima in the plot
Chiew2010 apap_sul looks good for the first time, even without scaling
load("methacetin_to_lung.rda")
x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)
# get the only the parameters needed for x
pars_to_lung <- rep(0,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
pars <- c(all_pars, pars_to_lung)[getParameters(x)]
free_parameters_to_lung <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km",
"Ka_apap",
paste0("k_to_lung_", c("sul", "", "cys", "glu"))
# "CLrenal_apap", "CLrenal_metc13", "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
)
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_to_lung <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_to_lung <- c(free_parameters_to_lung, scale_parameters_to_lung)
# error_model_to_lung <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung <- setdiff(getSymbols(error_model_to_lung), names(error_model_to_lung)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_to_lung] <- paste0("exp(log", free_parameters_to_lung, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung, x = scale_parameters_to_lung, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung, x = error_parameters_to_lung, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_to_lung <- NULL
for(i in 1:length(p_list)) { p_to_lung <<- p_to_lung + p_list[[i]]}
g_to_lung <- Y(observables_to_lung, x_to_lung)#, parameters = c(free_parameters_to_lung, scale_parameters_to_lung))
# err_to_lung <- Y(error_model_to_lung, g_to_lung)
obj_to_lung <- normL2(mydatalist, (g_to_lung*x_to_lung*p_to_lung))#, errmodel = err_to_lung)
pouter_to_lung <- rep(0, length(getParameters(obj_to_lung))) %>% set_names(getParameters(obj_to_lung))
pouter_to_lung[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
obj_to_lung(pouter_to_lung)
# job_to_lung <- runbg({myfit <- mstrust(objfun = obj_to_lung, center = pouter_to_lung, sd = 2, studyname = "methacetin", cores = 24, fits = 100, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_to_lung, file = "job_to_lung.rda")
# job_to_lung$check()# myfit_to_lung <- job_to_lung$get()$ruprecht1
# save(myfit_to_lung, file = "myfit_to_lung.rda")
# job_to_lung$purge()
load("myfit_to_lung.rda")
loadDLL(x_to_lung)## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.somyfit_to_lung %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung %>% as.parframe())+scale_y_log10()
mypred_to_lung <- (g_to_lung*x_to_lung*p_to_lung)(mytimes*4, myfit_to_lung %>% as.parframe() %>% as.parvec(1) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
mypred_to_lung2 <- (g_to_lung*x_to_lung*p_to_lung)(mytimes*4, myfit_to_lung %>% as.parframe() %>% as.parvec(55) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F) %>% set_names(paste(names(.),"_2"))
# myplot_to_lung <- plotCombined(mypred_to_lung, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung <- plotCombined(c(mypred_to_lung,mypred_to_lung2), mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung)
# myplot_to_lung## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 70 1557.189 TRUE 142 1.882061 5.800306
## 2 34 1557.199 TRUE 148 1.296540 5.215847
## 3 41 1557.213 TRUE 142 1.652053 5.570782
## 4 4 1557.214 TRUE 150 1.600278 5.520080
## 5 54 1557.249 FALSE 150 1.315475 5.237448
## 6 95 1557.256 FALSE 150 1.215461 5.132848
## logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 1.3139749 5.074520 -7.348043 -6.117347
## 2 2.0624390 5.822052 -7.348164 -6.117514
## 3 1.0046156 4.764574 -7.348081 -6.117433
## 4 1.0356908 4.794728 -7.348121 -6.117488
## 5 0.7347650 4.491662 -7.348263 -6.117822
## 6 0.7286464 4.489584 -7.348001 -6.117180
## logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1 -0.1459383 -7.460770 -0.2039797 3.060832
## 2 -0.1448628 -7.460689 -0.2039490 3.060762
## 3 -0.1454166 -7.460765 -0.2039929 3.060800
## 4 -0.1444364 -7.460705 -0.2039785 3.060779
## 5 -0.1423281 -7.460575 -0.2039709 3.060692
## 6 -0.1464887 -7.460877 -0.2040240 3.060828
## logk_to_lung_glu
## 1 -0.8190345
## 2 -0.8201550
## 3 -0.8195732
## 4 -0.8207004
## 5 -0.8229671
## 6 -0.8182854
## logAPAPCYS_HLM_CL: 0.00064385123805241
## logAPAPCYS_Km: 0.0022042952730362
## logAPAPGLU_HLM_CL: 6.56702557689903
## logAPAPGLU_Km: 330.40080787627
## logAPAPSUL_HLM_CL: 3.72093453867929
## logAPAPSUL_Km: 159.895379325816
## logKa_apap: 0.000575213236929382
## logk_to_lung_: 0.815478948433962
## logk_to_lung_cys: 21.3453185074998
## logk_to_lung_glu: 0.44085707897448
## logk_to_lung_sul: 0.864211035946343increase rate to lung and allow scaling - didn’t work out, many fits didn’t converge
The plots look OK, but the values don’t make any sense.
x_to_lung <- Xs(myodemodel_to_lung) # make prediction function
loadDLL(x_to_lung)
# get the only the parameters needed for x
pars_to_lung_2 <- rep(0,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
pars <- c(all_pars, pars_to_lung_2)[getParameters(x)]
free_parameters_to_lung_2 <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km",
"Ka_apap",
paste0("k_to_lung_", c("sul", "", "cys", "glu"))
# "CLrenal_apap", "CLrenal_metc13", "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
)
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_to_lung_2 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung_2 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_to_lung_2 <- c(free_parameters_to_lung_2, scale_parameters_to_lung_2)
# error_model_to_lung_2 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung_2 <- setdiff(getSymbols(error_model_to_lung_2), names(error_model_to_lung_2)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_to_lung_2] <- paste0("exp(log", free_parameters_to_lung_2, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_2, x = scale_parameters_to_lung_2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_2)
scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
# scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung_2, x = error_parameters_to_lung_2, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung_2)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_to_lung_2 <- NULL
for(i in 1:length(p_list)) { p_to_lung_2 <<- p_to_lung_2 + p_list[[i]]}
g_to_lung_2 <- Y(observables_to_lung_2, x)#, parameters = c(free_parameters_to_lung_2, scale_parameters_to_lung_2))
# err_to_lung_2 <- Y(error_model_to_lung_2, g_to_lung_2)
obj_to_lung_2 <- normL2(mydatalist, (g_to_lung_2*x_to_lung*p_to_lung_2))#, errmodel = err_to_lung_2)
pouter_to_lung_2 <- rep(0, length(getParameters(obj_to_lung_2))) %>% set_names(getParameters(obj_to_lung_2))
pouter_to_lung_2[names(myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec())] <- myfit5 %>% as.parframe() %>% {.[2,]} %>% as.parvec()
# obj_to_lung_2(pouter_to_lung_2)
# job_to_lung_2 <- runbg({myfit <- mstrust(objfun = obj_to_lung_2, center = pouter_to_lung_2, sd = 2, studyname = "methacetin", cores = 24, fits = 100, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung_2", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_to_lung_2, file = "job_to_lung_2.rda")
# job_to_lung_2$check()# myfit_to_lung_2 <- job_to_lung_2$get()$ruprecht1
# save(myfit_to_lung_2, file = "myfit_to_lung_2.rda")
# job_to_lung_2$purge()
load("myfit_to_lung_2.rda")
loadDLL(x_to_lung)## The following local files were dynamically loaded: methacetin_to_lung.so, methacetin_to_lung_s.somyfit_to_lung_2 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung_2 %>% as.parframe())+scale_y_log10()
mypred_to_lung_2 <- (g_to_lung_2*x_to_lung*p_to_lung_2)(mytimes*4, myfit_to_lung_2 %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_to_lung_2 <- plotCombined(mypred_to_lung_2, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung_2 <- plotCombined(mypred_to_lung_2, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung_2)
# myplot_to_lung_2## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 61 811.8280 FALSE 150 -6.998626 -17.50336
## 2 7 812.1984 FALSE 150 -7.019726 -16.97033
## 3 79 812.3622 FALSE 150 -7.031897 -16.94359
## 4 46 812.3956 FALSE 150 -7.033562 -16.92471
## 5 36 822.9495 FALSE 150 -7.920227 -17.47499
## 6 54 831.7163 FALSE 150 -9.112826 -16.08981
## logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 -4.647647 -7.067046 -8.005360 -25.35451
## 2 -4.677744 -7.095040 -8.017149 -24.29624
## 3 -4.699630 -7.149380 -8.015090 -24.25044
## 4 -4.703237 -7.161027 -8.014722 -24.18622
## 5 -4.908147 -8.078681 -7.579660 -23.58753
## 6 -5.018467 -8.363523 -7.356687 -23.10982
## logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1 -2.288579 -8.479607 -4.093847 -2.952564
## 2 -2.275876 -8.451360 -4.094141 -2.997499
## 3 -2.269219 -8.442444 -4.094653 -3.017023
## 4 -2.268317 -8.441194 -4.094752 -3.020282
## 5 -2.031370 -8.290519 -4.068860 -3.637812
## 6 -1.887041 -8.282598 -4.051058 -3.751590
## logk_to_lung_glu log_scale_apap_Chan1997_14_0
## 1 6.776023 2.035524
## 2 6.772015 1.963864
## 3 7.088785 1.967955
## 4 6.708936 1.971821
## 5 6.577030 2.122688
## 6 5.219719 2.144451
## log_scale_apap_glu_Chan1997_14_0 log_scale_apap_sul_Chan1997_14_0
## 1 6.495936 -2.067167
## 2 6.514004 -2.111232
## 3 6.844040 -2.125823
## 4 6.465985 -2.127919
## 5 7.263784 -2.351949
## 6 7.126559 -2.320297
## log_scale_apap_cys_Chan1997_14_0 log_scale_apap_Chiew2010_56_0
## 1 -4.984626 -0.1457100
## 2 -4.992699 -0.2872532
## 3 -5.003292 -0.3473295
## 4 -5.005029 -0.3556862
## 5 -5.625267 -1.0307136
## 6 -5.839268 -1.2249324
## log_scale_apap_glu_Chiew2010_56_0 log_scale_apap_sul_Chiew2010_56_0
## 1 8.738710 -2.787880
## 2 8.752954 -2.781700
## 3 9.079320 -2.776042
## 4 8.700771 -2.775228
## 5 9.406054 -2.621662
## 6 9.203844 -2.490344
## log_scale_apap_Critchley2005_14_0 log_scale_apap_glu_Critchley2005_14_0
## 1 1.914790 6.775920
## 2 1.843805 6.792860
## 3 1.848134 7.121887
## 4 1.852019 6.743747
## 5 2.008014 7.506618
## 6 2.029117 7.351299
## log_scale_apap_sul_Critchley2005_14_0
## 1 -2.517989
## 2 -2.561130
## 3 -2.575354
## 4 -2.577414
## 5 -2.786350
## 6 -2.748679
## log_scale_apap_cys_Critchley2005_14_0 log_scale_apap_Rawlins1977_0_1
## 1 -5.196222 -2.667753
## 2 -5.205343 -2.668026
## 3 -5.216401 -2.668411
## 4 -5.218228 -2.668487
## 5 -5.860658 -2.650692
## 6 -6.081074 -2.638401
## log_scale_apap_Rawlins1977_1_0 log_scale_apap_Rawlins1977_2_0
## 1 2.251798 1.204153
## 2 2.148307 1.137746
## 3 2.142583 1.140532
## 4 2.145176 1.143994
## 5 2.125534 1.249054
## 6 2.134132 1.174536
## log_scale_apap_Rawlins1977_05_0
## 1 6.104588
## 2 5.154957
## 3 4.939859
## 4 4.913994
## 5 3.014949
## 6 2.848154
## logAPAPCYS_HLM_CL: 0.000333669227703504
## logAPAPCYS_Km: 9.74264737617257e-12
## logAPAPGLU_HLM_CL: 0.00091313568705377
## logAPAPGLU_Km: 2.50258813575314e-08
## logAPAPSUL_HLM_CL: 0.00958412454265644
## logAPAPSUL_Km: 0.000852748222445177
## logKa_apap: 0.000207660238495109
## logk_to_lung_: 0.016674967733492
## logk_to_lung_cys: 0.0522056534562965
## logk_to_lung_glu: 876.576074026939
## logk_to_lung_sul: 0.101410478546347
## log_scale_apap_Chan1997_14_0: 7.65626070083065
## log_scale_apap_Chiew2010_56_0: 0.864408357893709
## log_scale_apap_Critchley2005_14_0: 6.78551630713505
## log_scale_apap_cys_Chan1997_14_0: 0.00684233640219463
## log_scale_apap_cys_Critchley2005_14_0: 0.00553744269024556
## log_scale_apap_glu_Chan1997_14_0: 662.443889437685
## log_scale_apap_glu_Chiew2010_56_0: 6239.84412853667
## log_scale_apap_glu_Critchley2005_14_0: 876.485712326581
## log_scale_apap_Rawlins1977_0_1: 0.0694080270944431
## log_scale_apap_Rawlins1977_05_0: 447.908262977161
## log_scale_apap_Rawlins1977_1_0: 9.50480772104851
## log_scale_apap_Rawlins1977_2_0: 3.33393329597214
## log_scale_apap_sul_Chan1997_14_0: 0.126543722411801
## log_scale_apap_sul_Chiew2010_56_0: 0.0615515731003974
## log_scale_apap_sul_Critchley2005_14_0: 0.0806215896552064Estimate dosing - some things look OK
Here I try to estimate the dosing. Dunno if this makes much sense, but it might be worth a try.
load("methacetin.rda")
x <- Xs(myodemodel) # make prediction function
loadDLL(x)
# get the only the parameters needed for x
pars <- all_pars[getParameters(x)]
free_parameters_est_input <- c("APAPGLU_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPSUL_HLM_CL", # Vmax value
"APAPGLU_Km", # Km value
"APAPCYS_HLM_CL", # Vmax value
"APAPCYS_Km", # Km value
"Ka_apap" ,
"D_apap",
"CLrenal_apap", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_apap_sul"
#, #"F_apap_sul" ,
# "Kpre_apap", "Kpki_apap", "Kpli_apap",
# "Kpre_apap_cys", "Kpki_apap_cys", "Kpli_apap_cys",
# "Kpre_apap_glu", "Kpki_apap_glu", "Kpli_apap_glu",
# "Kpre_apap_sul", "Kpre_apap_glu", "Kpli_apap_glu"#,
# "Kpre_co2c13", "Kpre_co2c13", "Kpli_co2c13",
# "Kpre_metc13", "Kpre_metc13", "Kpli_metc13"
)
fixed_parameters_est_input <- pars[!(names(pars)%in%c(free_parameters_est_input,names(f)[1]))] %>% names
mydatalist <- data %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap_cys") & time == 5400)) %>%
filter(!((study %>% str_detect("Chan")) & (name %in% "apap") & time == 18000)) %>%
fitErrorModel(factors = c("study", "D_apap", "Ave_apap", "name")) %>%
select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_est_input <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_est_input <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_est_input <- c(free_parameters_est_input, scale_parameters_est_input)
# error_model_est_input <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_est_input <- setdiff(getSymbols(error_model_est_input), names(error_model_est_input)) %>% set_names(.,.)
i <- 1
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_est_input] <- paste0("exp(log", free_parameters_est_input, ")")
if(cond[1]==0) {
trafo["Ave_apap"] <- paste0("exp(logAve_apap_", rownames(conditions)[i], ")") %>% str_replace_all("\\.", "")
trafo["D_apap"] <- "0"
} else {
trafo["D_apap"] <- paste0("exp(logD_apap_", rownames(conditions)[i], ")") %>% str_replace_all("\\.", "")
trafo["Ave_apap"] <- 0
}
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_2, x = scale_parameters_to_lung_2, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_2)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_est_input, x = error_parameters_est_input, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_est_input)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_est_input <- NULL
for(i in 1:length(p_list)) { p_est_input <<- p_est_input + p_list[[i]]}
g_est_input <- Y(observables_est_input, x)#, parameters = c(free_parameters_est_input, scale_parameters_est_input))
# err_est_input <- Y(error_model_est_input, g_est_input)
obj_est_input <- normL2(mydatalist, (g_est_input*x*p_est_input))#, errmodel = err_est_input)
pouter_est_input <- rep(0, length(getParameters(obj_est_input))) %>% set_names(getParameters(obj_est_input))
pouter_est_input[names(myfit_CLrenal %>% as.parframe() %>% as.parvec())] <- myfit_CLrenal %>% as.parframe() %>% as.parvec()
pouter_est_input <- pouter_est_input[order(names(pouter_est_input))]
pouter_est_input[c(6,11:16)] <- log(c(1,1.4,5.6,1.4,0.5,1,2))
# p_est_input(pouter_est_input)
# job_est_input <- runbg({myfit <- mstrust(objfun = obj_est_input, center = pouter_est_input, sd = 5, studyname = "methacetin", cores = 16, fits = 200); myfit}, machine = "ruprecht2", filename = "job_est_input", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_est_input, file = "job_est_input.rda")
# job_est_input$check()
# (g_est_input*x*p_est_input)(times = mytimes, pouter_est_input, deriv = F) %>% plotPrediction(name %in% names(observables))# load("job_est_input.rda")
# myfit_est_input <- job_est_input$get()$ruprecht2
# save(myfit_est_input, file = "myfit_est_input.rda")
# job_est_input$purge()
load("myfit_est_input.rda")
loadDLL(x)## The following local files were dynamically loaded: methacetin.so, methacetin_s.somyfit_est_input %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_est_input %>% as.parframe())+scale_y_log10()
mypred_est_input <- (g_est_input*x*p_est_input)(mytimes*4, myfit_est_input %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_est_input <- plotCombined(mypred_est_input, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_est_input <- plotCombined(mypred_est_input, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_est_input)
# myplot_est_input## index value converged iterations logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 58 1122.299 FALSE 100 -10.48928 -4.957525
## 2 156 1122.300 TRUE 90 -10.49027 -4.961741
## 3 91 1142.864 TRUE 54 -10.48757 -5.127342
## 4 192 1142.864 TRUE 77 -10.48709 -5.124875
## 5 130 1142.864 TRUE 43 -10.48786 -5.128740
## 6 144 1142.864 TRUE 75 -10.48767 -5.127869
## logAPAPGLU_HLM_CL logAPAPGLU_Km logAPAPSUL_HLM_CL
## 1 7.615842 10.65704 -5.418030
## 2 6.535269 9.57645 -5.418026
## 3 10.713608 13.62289 -5.292063
## 4 18.348451 21.25771 -5.292033
## 5 12.851771 15.76106 -5.292060
## 6 9.264502 12.17378 -5.292044
## logAve_apap_Rawlins1977_0_1 logCLrenal_apap logCLrenal_apap_cys
## 1 0.10147172 -6.716203 -5.779922
## 2 0.10146728 -6.716345 -5.778142
## 3 0.02988996 -34.131247 -5.648947
## 4 0.02989775 -32.298981 -5.650080
## 5 0.02989460 -31.947337 -5.648331
## 6 0.02989872 -28.053280 -5.648702
## logCLrenal_apap_glu logCLrenal_apap_sul logD_apap_Chan1997_14_0
## 1 -5.444322 -4.590098 0.8099026
## 2 -5.444325 -4.590114 0.8099012
## 3 -5.331795 -4.452475 0.7546452
## 4 -5.331763 -4.452424 0.7546769
## 5 -5.331793 -4.452482 0.7546558
## 6 -5.331771 -4.452459 0.7546687
## logD_apap_Chiew2010_56_0 logD_apap_Critchley2005_14_0
## 1 2.563208 0.6401644
## 2 2.563184 0.6401447
## 3 2.467657 0.5631522
## 4 2.467679 0.5631632
## 5 2.467656 0.5631472
## 6 2.467666 0.5631538
## logD_apap_Rawlins1977_05_0 logD_apap_Rawlins1977_1_0
## 1 -0.8721721 0.12574694
## 2 -0.8721339 0.12575422
## 3 -0.9409273 0.05866417
## 4 -0.9409427 0.05867302
## 5 -0.9409016 0.05867707
## 6 -0.9408988 0.05868636
## logD_apap_Rawlins1977_2_0 logKa_apap
## 1 0.8141646 -7.855334
## 2 0.8141561 -7.855325
## 3 0.7695203 -7.811703
## 4 0.7695401 -7.811767
## 5 0.7695220 -7.811715
## 6 0.7695283 -7.811742
## logAPAPCYS_HLM_CL: 2.78331452999271e-05
## logAPAPCYS_Km: 0.00703030850292055
## logAPAPGLU_HLM_CL: 2030.10329039433
## logAPAPGLU_Km: 42490.6189263482
## logAPAPSUL_HLM_CL: 0.00443587568913188
## logAve_apap_Rawlins1977_0_1: 1.1067986135081
## logCLrenal_apap: 0.00121112777686759
## logCLrenal_apap_cys: 0.0030889577568275
## logCLrenal_apap_glu: 0.00432076935688847
## logCLrenal_apap_sul: 0.0101518605119084
## logD_apap_Chan1997_14_0: 2.24768909406063
## logD_apap_Chiew2010_56_0: 12.9773760893221
## logD_apap_Critchley2005_14_0: 1.896792740859
## logD_apap_Rawlins1977_05_0: 0.418042515914904
## logD_apap_Rawlins1977_1_0: 1.13399516723761
## logD_apap_Rawlins1977_2_0: 2.25728919295756
## logKa_apap: 0.000387678645051781Do the estimation of _to_lung again with the original data instead of the data with fitted error model - Problem in many fits.
# Multiply rates from "ve" to "lu" with additional free parameters
reactions_to_lung_3 <- reactions
reactions_to_lung_3[reactions_to_lung_3$from %>% str_detect("(1*Ave_apap)"),3] %<>% paste0(" * (k_to_lung_", c("sul", "", "cys", "glu"),")")
reactions_to_lung_3
el_to_lung_3 <- eqnlist()
for(i in 1:nrow(reactions_to_lung_3)) el_to_lung_3 <- addReaction(el_to_lung_3, reactions_to_lung_3$from[i], reactions_to_lung_3$to[i], reactions_to_lung_3$rate[i], reactions_to_lung_3$description[i])
# Convert to "eqnvec", which is basically a named vector of the ODEs and the names denote the states
f_to_lung_3 <- el_to_lung_3 %>% as.eqnvec()
pars_to_lung_3 <- rep(1,4) %>% set_names(paste0("k_to_lung_", c("sul", "", "cys", "glu")))
parameters <- c(all_pars, pars_to_lung_3)[getParameters(x)] %>% names()
free_parameters_to_lung_3 <- parameters[parameters %>% sapply(. %>% str_detect(c("_CL", "Km", "Ka_apap$", "k_to")) %>% any) ]
# myodemodel_to_lung_3 <- odemodel(f_to_lung_3, modelname = "methacetin_to_lung_3", fixed = setdiff(parameters, c(free_parameters_to_lung_3, "Ave_apap"))) # Include "Ave_apap" because of a bug: You need the sensitivities to at least one initial value
# save(myodemodel_to_lung_3, file = "methacetin_to_lung_3.rda")
load("methacetin_to_lung_3.rda")
x_to_lung_3 <- Xs(myodemodel_to_lung_3) # make prediction function
loadDLL(x_to_lung_3)
# get the only the parameters needed for x
pars <- c(all_pars, pars_to_lung_3)[getParameters(x)]
free_parameters_to_lung_3 <- c("APAPGLU_HLM_CL", "APAPGLU_Km", "APAPSUL_HLM_CL", "APAPSUL_Km", "APAPCYS_HLM_CL", "APAPCYS_Km",
"Ka_apap",
paste0("k_to_lung_", c("sul", "", "cys", "glu"))
# "CLrenal_apap", "CLrenal_metc13", "CLrenal_apap_sul", "CLrenal_apap_cys", "CLrenal_apap_glu", "CLrenal_co2c13"
)
mydatalist <- data %>% filter(!is.na(sigma)) %>% select(-n) %>% as.datalist()
conditions <- mydatalist %>% attr("condition.grid")
observables_to_lung_3 <- c(apap = "Ave_apap/(BW*FVve)*scale_apap",
apap_glu = "Ave_apap_glu/(BW*FVve)*scale_apap_glu",
apap_sul = "Ave_apap_sul/(BW*FVve)*scale_apap_sul",
apap_cys = "Ave_apap_cys/(BW*FVve)*scale_apap_cys")
scale_parameters_to_lung_3 <- paste0("scale_apap", c("", "_glu", "_sul", "_cys")) %>% set_names(.,.)
# free_parameters_to_lung_3 <- c(free_parameters_to_lung_3, scale_parameters_to_lung_3)
# error_model_to_lung_3 <- c(apap = "srel_apap*apap^2 +s0_apap",
# apap_glu = "srel_apap_glu*apap_glu^2 +s0_apap_glu",
# apap_sul = "srel_apap_sul*apap_sul^2 +s0_apap_sul",
# apap_cys = "srel_apap_cys*apap_cys^2 +s0_apap_cys")
# error_parameters_to_lung_3 <- setdiff(getSymbols(error_model_to_lung_3), names(error_model_to_lung_3)) %>% set_names(.,.)
p_list <- lapply(1:nrow(conditions), function(i) {
cond <- unlist(conditions[i,])[2:3]
trafo <- as.character(pars) %>% set_names(names(pars))
trafo[names(cond)] <- cond
trafo[free_parameters_to_lung_3] <- paste0("exp(log", free_parameters_to_lung_3, ")")
scales <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", scale_parameters_to_lung_3, x = scale_parameters_to_lung_3, y = .)} %>% str_replace_all("\\.", "") %>% set_names(scale_parameters_to_lung_3)
# scales <- scales[names(scales) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
scales[1:length(scales)] <- "1"
# errors <- rownames(conditions)[i] %>% {repar("x~exp(log_x_y)", error_parameters_to_lung_3, x = error_parameters_to_lung_3, y = .)}%>% str_replace_all("\\.", "") %>% set_names(error_parameters_to_lung_3)
# errors <- errors[names(errors) %>% sapply(. %>% str_detect(mydatalist[[i]][["name"]] %>% unique() %>% paste0("$")) %>% any)]
trafo <- c(trafo, scales)#, errors)
p <- P(trafo, condition=rownames(conditions[i,]))
return(p)
})
p_to_lung_3 <- NULL
for(i in 1:length(p_list)) { p_to_lung_3 <<- p_to_lung_3 + p_list[[i]]}
g_to_lung_3 <- Y(observables_to_lung_3, x_to_lung_3)#, parameters = c(free_parameters_to_lung_3, scale_parameters_to_lung_3))
# err_to_lung_3 <- Y(error_model_to_lung_3, g_to_lung_3)
obj_to_lung_3 <- normL2(mydatalist, (g_to_lung_3*x_to_lung_3*p_to_lung_3))#, errmodel = err_to_lung_3)
getParameters(obj_to_lung_3)
myfit_to_lung %>% as.parframe() %>% as.parvec()
pouter_to_lung_3 <- myfit_to_lung %>% as.parframe() %>% as.parvec()
# job_to_lung_3 <- runbg({myfit <- mstrust(objfun = obj_to_lung_3, center = pouter_to_lung_3, sd = 4, studyname = "methacetin", cores = 24, fits = 1000, iterlim = 150); myfit}, machine = "ruprecht1", filename = "job_to_lung_3", input = global_env_without(c("mypred", "job", "myfit", "myplot")))
# save(job_to_lung_3, file = "job_to_lung_3.rda")
# job_to_lung_3$check()# load("job_to_lung_3.rda")
# myfit_to_lung_3 <- job_to_lung_3$get()$ruprecht1
# save(myfit_to_lung_3, file = "myfit_to_lung_3.rda")
# job_to_lung_3$purge()
load("myfit_to_lung_3.rda")
loadDLL(x_to_lung_3)## The following local files were dynamically loaded: methacetin_to_lung_3.so, methacetin_to_lung_3_s.somyfit_to_lung_3 %>% as.parframe() %T>% {head(.) %>% print} %>% as.parvec() %>% exp()
plotValues(myfit_to_lung_3 %>% as.parframe())+scale_y_log10()
mypred_to_lung_3 <- (g_to_lung_3*x_to_lung_3*p_to_lung_3)(mytimes*4, myfit_to_lung_3 %>% as.parframe() %>% as.parvec(35) %>% {names(.) <- names(.) %>% str_replace_all("\\.","");.}, deriv = F)
# myplot_to_lung_3 <- plotCombined(mypred_to_lung_3, mydatalist, name %>% sapply(. %>% str_detect(c(paste0("^",names(observables)), "Agu", "Ali", "Ave")) %>% any))
myplot_to_lung_3 <- plotCombined(mypred_to_lung_3, mydatalist, name %in% names(observables))
plotly::ggplotly(myplot_to_lung_3)
# myplot_to_lung_3## index value converged iterations logAPAPGLU_HLM_CL logAPAPGLU_Km
## 1 598 1502.666 FALSE 150 5.222239 9.111693
## 2 585 1502.666 FALSE 150 6.010086 9.899766
## 3 230 1502.691 FALSE 150 7.145232 11.034067
## 4 74 1502.694 FALSE 150 3.539114 7.429065
## 5 901 1502.698 FALSE 150 3.013285 6.905417
## 6 852 1502.701 FALSE 150 2.820840 6.710850
## logAPAPSUL_HLM_CL logAPAPSUL_Km logAPAPCYS_HLM_CL logAPAPCYS_Km
## 1 4.908114 9.176861 -7.698137 -6.722977
## 2 3.700609 7.968713 -7.698307 -6.723640
## 3 1.294309 5.563838 -7.698262 -6.723571
## 4 2.076930 6.344381 -7.698409 -6.723928
## 5 3.797133 8.063931 -7.696645 -6.715827
## 6 5.819994 10.087409 -7.698370 -6.723744
## logk_to_lung_sul logKa_apap logk_to_lung_ logk_to_lung_cys
## 1 -0.4319789 -7.093422 -0.05287355 2.874983
## 2 -0.4313477 -7.093479 -0.05289277 2.874915
## 3 -0.4328742 -7.093376 -0.05285579 2.874963
## 4 -0.4307210 -7.093524 -0.05290436 2.874860
## 5 -0.4297375 -7.093220 -0.05277499 2.875319
## 6 -0.4306401 -7.093529 -0.05290571 2.874868
## logk_to_lung_glu
## 1 -0.7626255
## 2 -0.7628978
## 3 -0.7619719
## 4 -0.7632058
## 5 -0.7651345
## 6 -0.7632765
## logAPAPCYS_HLM_CL: 0.000453671805412169
## logAPAPCYS_Km: 0.00120295170864251
## logAPAPGLU_HLM_CL: 185.348656922024
## logAPAPGLU_Km: 9060.6175670097
## logAPAPSUL_HLM_CL: 135.383862495601
## logAPAPSUL_Km: 9670.74882908226
## logKa_apap: 0.000830550499923049
## logk_to_lung_: 0.948499944225183
## logk_to_lung_cys: 17.7251183888039
## logk_to_lung_glu: 0.466440188442027
## logk_to_lung_sul: 0.649223078801294# save(list= ls(), file = "workspace.rda")Next possibilities
- Take out some data, eg. apap to fit the other parameters such that the fit isn’t dominated by apap (not so good)
- Estimate the Km and Vmaxes each individually, then use those values as a prior? (maybe that’s not so good, because: what are the values of the other Kms and Vmaxes?)